System modeling, control and optimization

ABSTRACT

A method for controlling an operation of a system that includes a disturbance rejection model that is configured for modeling the operation of the system so to generate a predicted value for a system output at a future time. The disturbance rejection model may include a neural network for mapping system inputs to the system output. The method may include the steps of: calculating a probabilistic distribution for the predicted value of the system output at the future time, where the calculating of the probabilistic distribution includes a Bayesian evidence framework without sampling; and controlling the operation of the system per the calculated probabilistic distribution.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to U.S. Provisional Patent ApplicationNo. 62/366,248 entitled “PROCESS MODELING AND CONTROL USING MODELS WITHCALCULATED CONFIDENCE” filed on Jul. 25, 2016; which provisionalapplication is incorporated herein by reference in its entirety; thisapplication claims the benefit of the provisional's filing date under 35U.S.C. 119(e).

BACKGROUND OF THE INVENTION

The present invention relates to systems and methods of process orsystem modeling and control and, specifically, process modeling andcontrol using a calculated confidence in the system model.

In a conventional fossil fuel-fired power plant or generating unit, afossil fuel/air mixture is ignited in a boiler. Large volumes of waterare pumped through tubes inside the boiler, and the intense heat fromthe burning fuel turns the water in the boiler tubes into high-pressuresteam. In an electric power generating application, the high-pressuresteam from the boiler passes into a turbine comprised of a plurality ofturbine blades. Once the steam hits the turbine blades, it causes theturbine to spin rapidly. The shaft of the spinning turbine is linked tothe shaft of a generator, and the rotating shaft within the generatormay be used to create an electric potential.

As used herein, the term “power generating plant” refers to one or morepower generating units. Each power generating unit may drive one or moreturbines used for generating electricity. A power generating unit istypically powered by fossil fuels (including but not limited to, coal,natural gas or oil), and includes a boiler for producing hightemperature steam; air pollution control (APC) devices for removal ofpollutants from flue gas; a stack for release of flue gas; and a watercooling system for condensing the high temperature steam. A typicalpower generating unit will be described in detail below in relation toFIG. 1.

As will be appreciated, boiler combustion or other characteristics of afossil fuel-fired power generating unit are influenced by dynamicallyvarying parameters of the power generating unit, including, but notlimited to, air-to-fuel ratios, operating conditions, boilerconfiguration, slag/soot deposits, load profile, fuel quality andambient conditions. Changes to the business and regulatory environmentshave increased the importance of dynamic factors such as fuelvariations, performance criteria, emissions control, operatingflexibility and market driven objectives (e.g., fuel prices, cost ofemissions credits, cost of electricity, etc.).

Further, about one half of the electric power generated in the UnitedStates is generated using coal-fired power generating units. Coal-firedpower generating units used in power plants typically have an assortmentof air pollution control (APC) devices installed for reducing nitrogenoxides (NOx), sulfur oxides (SOx), and particulate emissions. In thisregard, selective catalytic reduction (SCR) systems are used for NOxreductions. Spray dry absorbers (SDA) and wet flue gas desulfurization(FGD) systems are used for SOx reductions. Electro-static precipitators(ESPs) and fabric filters (FF) are used for reducing particulateemissions.

Over the past two decades, combustion optimization systems have beenimplemented for advanced control of the combustion process within thefurnace. Typically, combustion optimization systems interface with thedistributed control system (DCS) of a power generating unit. Based uponthe current operating conditions of the power generating unit, as wellas a set of operator specified goals and constraints, the combustionoptimization system is used to compute the optimal fuel-to-air stagingwithin the furnace to achieve the desire goals and constraints.

Combustion optimization systems were originally implemented to reducenitrogen oxides (NOx) produced in the furnace and emitted to theatmosphere via the stack. For example, U.S. Pat. No. 5,280,756 teaches amethod and system for controlling and providing guidance in reducing NOxemissions based upon controllable combustion parameters and modelcalculations while maintaining satisfactory plant performance. U.S. Pat.No. 5,386,373 teaches the use of a predictive model of emissions,including NOx, in conjunction with a control system, while U.S. Pat. No.6,381,504 describes a method for optimally determining the distributionof air and fuel within a boiler by aggregating the distributions of airand fuel into two common variables, performing an optimization, and thencomputing the optimal distribution of fuel and air based upon theoptimal values of the aggregated variables. U.S. Pat. No. 6,712,604describes a system for controlling the combustion of fuel and air in aboiler such that the distributions of NOx and CO are maintained toaverage less than the maximum permitted levels.

In addition, combustion optimization approaches have been used tocontrol boiler parameters in addition to NOx, including unit heat rate,boiler efficiency, and mercury emissions. For example, U.S. Pat. No.7,398,652 teaches an approach to modeling controllable losses in a powergenerating unit and a method for optimizing the combustion process basedupon these controllable losses. U.S. Pat. No. 8,644,961 teaches a methodfor reducing mercury emissions from a coal-fired power plant whileobserving limits on the amount of carbon in the fly ash produced. U.S.Pat. No. 7,522,963 teaches a controller for directing operation of anair pollution control system, such as an FGD or SCR, such that apredefined optimization objective is minimized. Optimization techniqueshave also been used to control the removal of soot within a boiler. Forexample, U.S. Pat. No. 6,736,089 teaches a method for removal of sootbased upon optimizing a set of boiler performance parameters using amodel of the cleanliness factors of heat transfer surfaces. Optimizationhas also been extended beyond specific components within a powergenerating unit, such as the boiler, FGD and SCR. As provided in U.S.Pat. No. 7,844,351, which is hereby incorporated herein in its entirety,an approach to optimization of multiple components, within a singlepower generating unit or multiple power generating units, is described.

Neural networks have been used to both predict and optimize industrialsystems. For example, U.S. Pat. No. 5,167,009 teaches a method forpredicting output data of a process using a neural network andsubsequently using the prediction in a control or optimization system.As another example, U.S. Pat. No. 7,123,971 teaches an approach totraining neural network models to predict the change in the output basedupon the change in the inputs. Using this approach, referred to asdisturbance rejection based training, the model is trained to learn thecause of the change in an output rather than simply learning thecorrelation structure of the data which results in a significantimprovement in optimization results. In addition, U.S. Pat. No.6,725,208 teaches the use of a Bayesian training technique of combiningmultiple weighted neural network models to form a model used inoptimization of a process. This type of model can be used to estimatethe uncertainty in the predicted output, though these capabilities arelimited.

Thus, the prior art describes methods and systems for using models,including neural network models, to optimize industrial systems andprocesses. In addition, disturbance rejection based training has beenfound to provide more accurate neural network models for use in processoptimization. The prior art, however, does not describe a disturbancerejection based neural network training method for predicting change inthe output due to change in the inputs along with an estimate of theuncertainty in the predicted change in the output nor does it describehow such a neural network model could be used to improve optimization.The present invention provides systems and methods that overcome theabovementioned shortcomings of the prior art, and provide advantagesover conventional approaches to control and optimization industries suchas power generation.

BRIEF DESCRIPTION OF THE INVENTION

The present application thus describes a method for controlling anoperation of a system that includes a disturbance rejection model thatis configured for modeling the operation of the system so to generate apredicted value for a system output at a future time. The disturbancerejection model may include a neural network for mapping system inputsto the system output. The method may include the steps of: calculating aprobabilistic distribution for the predicted value of the system outputat the future time, where the calculating of the probabilisticdistribution includes a Bayesian evidence framework without sampling;and controlling the operation of the system per the calculatedprobabilistic distribution.

These and other features of the present application will become apparentupon review of the following detailed description of the preferredembodiments when taken in conjunction with the drawings and the appendedclaims.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features of this invention will be more completelyunderstood and appreciated by careful study of the following moredetailed description of exemplary embodiments of the invention taken inconjunction with the accompanying drawings, in which:

FIG. 1 shows a schematic representation of the primary components of aconventional fossil fuel power generating unit.

FIG. 2 illustrates a block diagram of an optimization system havingaspects in accordance with exemplary embodiments of the presentinvention.

FIG. 3 shows an example of the results of a “Design of Experiments” forcollecting data for training a model.

FIG. 4 illustrates an exemplary schematic diagram of an empirical model.

FIG. 5 shows a disturbance rejection model in accordance with aspects ofexemplary embodiments of the present invention.

FIG. 6 shows a flow diagram of the Bayesian disturbance rejectionalgorithm as may be used to compute the trained values of thehyperparameters, the weights, and the feedback coefficient of the modelin FIG. 5 in accordance with aspects of the present invention.

FIG. 7 shows the probability sample for a given point in the dataset asmay be computed pursuant to aspects of the present invention.

FIG. 8 illustrates how model confidence may be to expand the collectionof data used to train a model in accordance with embodiments of thepresent invention.

FIG. 9 shows a table comparing the results of the confidence calculationmay be compared to a threshold for determining whither additional dataneeds collection in accordance with embodiments of the presentinvention.

FIG. 10 illustrates a flow diagram according to the present inventionfor updating or training a model where additional data is collectedperiodically.

FIG. 11 illustrates a form of the disturbance rejection model that maybe used for optimization.

DETAILED DESCRIPTION OF THE INVENTION

The main components of a typical fossil fuel power generating unit 200will now be described with reference to FIG. 1. As illustrated, powergenerating unit 200 includes one or more forced draft (FD) fans 210 thatare powered by motors M1. Forced draft fans 210 supply air to mills 214and to burners 222, via an air preheater 212. Ambient air is heated asit passes through air preheater 212. Mills 214 include pulverizers thatare powered by motors M2. The pulverizers grind coal (or other fuel)into small particles (i.e., powder). The air received by the mills fromforced draft fans 210 is used to dry and carry the coal particles toburners 222. Air from forced draft fans 210 that is supplied to burners222, via air preheater 212, facilitates combustion of the coal atfurnace 224. Hot flue gas is drawn out of furnace 224 by one or moreinduced draft (ID) fans 260, and delivered to the atmosphere though achimney or stack 290. Induced draft fans 260 are powered by motors M3.Water is supplied to a drum 226 by control of a feedwater valve 228. Thewater in drum 226 is heated by furnace 224 to produce steam. This steamis further heated in a superheat region 230 by a superheater (notshown). A superheater spray unit (not shown) can introduce a smallamount of water to control the temperature of the superheated steam. Atemperature sensor (not shown) provides a signal indicative of thesensed temperature of the superheated steam. Superheated steam producedby power generating unit 200 is supplied to a turbine 250 that is usedto produce electricity. Steam received by the turbine is reused bycirculating the steam through a reheater (not shown) that reheats thesteam in a reheat region 240. A reheater spray unit (not shown) canintroduce a small amount of water to control the temperature of thereheated steam. A temperature sensor (not shown) provides a signalindicative of the sensed temperature of the reheated steam.

A “boiler” includes, but is not limited to, burners 222, furnace 224,drum 226, superheater, superheater spray unit, reheater, reheater sprayunit, mills 214, and a boiler economizer (not shown). The boilereconomizer recovers “waste heat” from the boiler's hot stack gas andtransfers this heat to the boiler's feedwater.

Soot cleaning devices (not shown), include, but are not limited to,sootblowers, water lances, and water cannons or hydro-jets. Sootcleaning devices use steam, water or air to dislodge deposits, such asslag, and clean surfaces throughout various locations in the boiler.Soot cleaning is required to maintain performance and efficiency ofpower generating unit 200. The number of soot cleaning devices on agiven power generating unit can range from several to over a hundred.Furthermore, the soot cleaning devices may be grouped together bylocation (e.g., zones in the boiler). Each group of soot cleaningdevices may be comprised of one or more soot cleaning devices. Forexample, a boiler may have eight soot cleaning device groups, each groupincluding five individual soot cleaning devices.

In addition, power generating unit 200 includes some form ofpost-combustion air pollution control (APC) equipment for removingpollutants from the flue gas. The APC equipment may include, but is notlimited to, a selective catalytic reactor (SCR) 206, an electro-staticprecipitator (ESP) 270, a fabric filter (FF) 272, a spray dry absorber(SDA) 274, and a wet flue gas desulfurization (FGD) system 276.

A selective catalytic reactor (SCR) is used to remove nitrogen oxides(NOx) from the flue gas. Dirty flue gas leaves the boiler and enters theselective catalytic reduction (SCR) system. Prior to entering the SCR,NOx in the inlet flue gas is measured with one or more analyzers. Inaddition, prior to entering the SCR, the flue gas passes through anammonia (NH₃) injection grid (not shown) located in the ductwork.Ammonia that has been mixed with dilution air is dosed into the flue gasby the injection grid. The NOx laden flue gas, ammonia and dilution airpass into the SCR reactor and over the SCR catalyst. The SCR catalystpromotes the reduction of NOx with ammonia to nitrogen and water. NOx“free” flue gas leaves the SCR reactor and exits the power generatingunit via potentially other APC subsystems and the stack.

Additional NOx analyzers are located in the NOx “free” flue gas streamexiting the SCR system or in the stack. The measured NOx outlet valueand the measured NOx inlet value are used to calculate a NOx removalefficiency. NOx removal efficiency is defined as the percentage of inletNOx removed from the flue gas.

In addition, a small amount of unreacted ammonia (i.e., “ammonia slip”)is exhausted from the SCR. This ammonia slip can react with othercomponents of the flue gas to form salts that can be deposited, andsubsequently foul other system components, such as the air preheater.Thus, to prevent fouling of components, the level of ammonia slip isoften constrained.

As the amount of ammonia injected into the flue gas increases, theremoval efficiency improves while the ammonia slip increases. Thus, aconstraint on ammonia slip indirectly constrains the removal efficiencyof the SCR. Because ammonia slip is often not directly measured on-linein real-time, it is typically indirectly controlled by limiting theremoval efficiency of the SCR.

An electro-static precipitator (ESP) is the most common approach toremoval of particulate matter from the flue gas steam of a powergenerating unit. In an ESP, particles suspended in the flue gas areelectrically charged. An electric field then forces the chargedparticles to an electrode where they are collected. A rapping system isused to remove the particles from the electrode. The removed particlesfall into an ash handle system which is used to dispose of the ash.Using this approach, ESPs can typically achieve 90%-99.5% removal ratesof particulate matter.

An ESP is typically comprised of a series of electrical plates withwires between the plates. The wires are used to charge the particlesusing corona discharge. An electric field for driving the particles isestablished between the wire and plates. The flue gas flows through aseries of electrically separated fields of plates and wires. Each ofthese fields may be separately powered. The primary motivation for usingseparate fields is to provide redundancy in the system.

A wet flue gas desulfurization (FGD) is the most common approach toremoval of significant amounts of SO₂ from the flue gas of powergenerating units. In a power generating unit, dirty, SO₂ laden flue gasis exhausted from a boiler. The SO₂ laden flue gas is input into anabsorber tower, which is the primary component in an FGD.

The SO₂ in the flue gas has a high acid concentration. Accordingly, theabsorber tower operates to place the SO₂ laden flue gas in contact witha liquid slurry having a higher pH level than that of the flue gas. Thisis accomplished by spraying the liquid slurry in countercurrent to theflue gas in the absorber tower.

During processing in the countercurrent absorber tower, the SO₂ in theflue gas will react with a calcium carbonate-rich slurry (limestone andwater) to form calcium sulfite, which is basically a salt and therebyremoving the SO₂ from the flue gas. The spray, including the SO₂ in theform of calcium sulfite, falls into a large tank at the bottom of theabsorber. The SO₂-cleaned flue gas is exhausted from the absorber tower,either to an exhaust stack or to downstream processing equipment. Ablower pressurizes ambient air to create oxidation air within theabsorber tank. The oxidation air is mixed with the slurry in the tank tooxidize the calcium sulfite to calcium sulfate. Each molecule of calciumsulfate binds with two molecules of water to form a compound that iscommonly referred to as gypsum. The gypsum is removed from the wet FGDprocessing unit and sold to, for example, manufacturers of constructiongrade wallboard. In order to sell the gypsum, it must be of anacceptable purity. The purity is affected by the pH which also affectsthe removal efficiency.

In FIG. 1, coal is used as the fuel for power generating unit 200. Ingeneral, fossil fuel power generating units often use a blend ofmultiple fuels. For example, most operators of coal-fired powergenerating units combine various types of coals to achieve a desiredblend that is burned in the furnace. Typically, several different typesof coal are stocked in the coal yard at a power generating plant. Thesedifferent coals may come from the same mine or from a variety of mines.If these coals are from the same mine, they may come from differentseams or different locations in the mine. Thus, each of the coals at thepower generating plant may have different costs, availability, and coalcharacteristics including heat, sulfur, nitrogen and ash content.Typically, the different coals are blended together by an operator oftenusing “rules of thumb” to supply the furnace with a desired blend ofcoal. In addition, fuel additives may be introduced into the blend toimprove heat rate or provide desired fuel characteristics.

It should be apparent that a typical power generating unit also includesadditional components well known to those skilled in the art, including,but not limited to, tubes for carrying fluids, valves, dampers, windbox,sensing devices for sensing a wide variety of system parameters (e.g.,temperature, pressure, flow rate, and flue gas components), andactuators for actuating components such as valves and dampers.

FIG. 2 illustrates a block diagram of an optimization system 100. In theillustrated embodiment, optimization system 100 is comprised of anoptimizer 110 and a model 120. Optimizer 110 and model 120 are bothdescribed in greater detail below. In accordance with an illustratedembodiment, optimization system 100 may form part of a supervisorycontroller 160 that communicates with a DCS 150. DCS 150 may be acomputer-based control system that provides regulatory control of apower generating plant 170. DCS 150 may take the form of a programmablelogic controller (PLC). Supervisory controller 160 is a computer systemthat provides supervisory control data to DCS 150. It should beunderstood that in an alternative embodiment, model 120 may reside on adifferent computer system than optimizer 110. An operator interface (notshown) may provide means for an operator to communicate with DCS 150.DCS 150 may also communicate with a historian (not shown).

Plant 170 may include one or more power generating units 200. Each powergenerating unit 200 may include a plurality of actuators 205 and sensors215. Actuators 205 may include devices for actuating components such asvalves and dampers. Sensors 215 may include devices for sensing varioussystem parameters (e.g., temperature, pressure, flow rate, and flue gascomponents).

As will be appreciated, model 120 is used to represent the relationshipbetween (a) manipulated variables (also “MV”) and disturbance variables(also “DV”) and (b) controlled variables (also “CV”). Manipulatedvariables (MVs) represent those variables that may be changed by theoperator or optimization system 100 to affect the controlled variables.The MVs typically include a subset of the following variables:over-fired air damper position, tilt and yaw biases, secondary airdamper position biases, fuel mill speed biases, primary air biases, O2bias, master burner tilt, windbox to furnace bias and other variables inthe boiler that may affect the CVs. As used herein, disturbancevariables refer to variables (associated with power generating unit 200)that affect the controlled variables, but cannot be manipulated orcontrolled by an operator (e.g., ambient conditions, characteristics ofthe coal, load, etc.). As will be appreciated, optimizer 110 mayfunction by determining an optimal set of setpoint values for themanipulated variables given (1) a desired goal associated with operationof the power generating unit (e.g., minimizing NOx production) and (2)constraints associated with operation of the power generating unit(e.g., limits on emissions of NOx, SO₂, CO₂, CO, mercury, ammonia slipand particulate matter).

At a predetermined frequency, optimization system 100 may obtain thecurrent values of manipulated variables, controlled variables anddisturbance variables from DCS 150. An “optimization cycle” commenceseach time the current values for the manipulated variables, controlledvariables and disturbance variables are read out from DCS 150. As willbe described in further detail below, optimization system 100 may usemodel 120 to determine an optimal set of setpoint values for themanipulated variables based upon current conditions of power generatingunit 200. The optimal set of setpoint values may then be sent to DCS150. An operator of plant 170 has the option of using the optimal set ofsetpoint values for the manipulated variables. In most cases, theoperator allows the computed optimal set of setpoint values for themanipulated variables to be used as setpoints values for control loops.Optimization system 100 runs in a closed loop adjusting the setpointsvalues of the manipulated variables at a predetermined frequency (e.g.,as frequently as every 10 seconds or as infrequently as every half hour)depending upon current operating conditions of power generating unit200.

Probabilistic Neural Network with Disturbance Rejection

The model of 120 may be developed based upon 1) known first principleequations describing the system, 2) data, resulting in an empiricalmodel, or 3) a combination of known first principle equations and data.In developing models for NOx and CO in a power plant, the firstprinciples equations are not easily derived—computational fluid dynamicsmust be used. For this reason, it is preferred to build empirical modelbased upon data collected from the plant. Furthermore, since the MVs mayhave been kept constant in the past or moved to control the unit, it istypical necessary to perform a “Design of Experiments” to collect thedata required to build an empirical model.

FIG. 3 shows an example of the results of a Design of Experiments forcollecting data for training a model. As shown, a time-series plot oftwo MVs, a DV, and a CV of interest (for example, NOx) are shown. MV₁,for example, may represent an upper separated over-fired air damper biasposition while MV₂, for example, may represent a lower separatedover-fired air damper bias position. The DV, for example, may representthe mega-watt load of the generator. In this illustration, it will beappreciated that, while only 2 MVs, 1 DV and 1 CV are shown, data istypically collected for many MVs, DVs, and CVs when performing a Designof Experiments collection of data. As further shown, the sampling timeof the data is shown at the top of the plot. In this case, samplingstarts at 11:00 am on Day 1 and concludes at 11:00 am on Day 2, withsamples being collected every fifteen minutes. The first sample islabeled time t=0, the second is labeled time t=1, etc., up to the lastpoint of time t=T.

To perform an appropriate Design of Experiments, the MVs may be movedindependently such that the moves are uncorrelated over time. It will beappreciated that the goal of the Design of Experiments is to collectsufficient data to develop the model, which then can be used foroptimization. As shown in FIG. 3, at each sampling time, the MVs aremoved (if they are to be moved for the sample period). After the sampletime, the MVs are held constant until the next sampling time. In thiscase, the interest is in building a steady state model and, thus, thesampling time is sufficient to allow the process to come to steadystate. Thus, for example, because a move in a damper position bias takesapproximately 15 minutes for its full effect to be realized, thesampling time is set accordingly. In the present application, the focusis primarily on steady state models. It should be appreciated, though,that present systems and methods may be extended to dynamic models. Asshown in FIG. 3, the effects of moving MV₁ and holding it steady fromtime (t=1) to just prior to time (t=2) influences the NOx at time (t=2).

The sampled points shown in FIG. 3 of the MVs, DV and NOx then may beused to create a dataset for training an empirical model, such as themodel shown in FIG. 4. As illustrated, the model in FIG. 4 (representedby the box M) has three inputs: MV₁(t−1), MV₂(t−1) and DV(t−1). As willbe appreciated, these inputs are used to predict the output NOx(t). Thatis, the inputs, the two MVs and DVs at time (t−1) are used to predictthe NOx at time (t), NO_(x) ^(P)(t). To properly capture therelationship between the manipulated/disturbance variables and thecontrolled variables, the model in FIG. 4 may need to be nonlinear. Thisis because nonlinear models can represent curved rather thanstraight-line relationships between manipulated/disturbance variablesand controlled variables, which are common to complex systems such asthose discussed herein. For example, a nonlinear, curved relationship istypically observed between over-fire air dampers and NOx levels. Giventhe foregoing requirements, a neural network based approach is presentlya preferred embodiment for implementing models in accordance with thepresent invention. Neural networks are developed based upon empiricaldata using advanced regression algorithms. See, for example, C. Bishop,Pattern Recognition and Machine Learning, Springer, New York, N.Y.,2006, fully incorporated herein by reference. Neural networks arecapable of capturing the nonlinearity commonly exhibited by boilers andother power generating systems.

As will be further appreciated, the model of FIG. 4 is independent ofthe current reading of NOx at time (t−1). Hence, it does not takeadvantage of the most current reading of NOx at time (t−1) to predictthe NOx at time (t). As described in U.S. Pat. No. 7,123,971, which isincorporated herein by reference in its entirety, using the value of NOxat time (t−1) may be used to significantly improve the accuracy of themodel.

FIG. 5 shows a disturbance rejection type model, which was described inU.S. Pat. No. 7,123,971. As shown, a neural network model (NN) is usedto predict the NOx at time (t−1) using the MVs and DVs at time (t−2).The difference between the prediction of the NOx at time (t−1) and theactual value of the NOx at time (t−1) are multiplied by a feedbackcoefficient (k) and then added to the prediction of NOx from the neuralnetwork at time (t). If the feedback coefficient is about 1, it will beappreciated that the prediction at time (t) is biased by the error inthe prediction at time (t−1).

The model of FIG. 5 has shown usefulness for optimization processes. Inoptimization, the interest is in minimizing the NOx at time (t+1) bydetermining a set of values for the MVs at time (t). The error in themodel at the current time is known and can be used to bias theprediction at time (t+1). As described in detail in U.S. Pat. No.7,123,971, this approach allows the rejection of slowly varyingunmeasured disturbance in the model prediction. Using the trainingmethod given in U.S. Pat. No. 7,123,971, more accurate models foroptimization can be determined. For this reason, the model of FIG. 5 isa preferred neural network embodiment.

U.S. Pat. No. 6,725,208, which is incorporated herein by reference inits entirety, introduces the use of a probabilistic model for use in anoptimization system. A probabilistic model may be provided that is basedupon Bayesian techniques for the model “M” shown in FIG. 4. Thisprobabilistic model is based upon a Bayesian sampling approach and,thus, uses a committee of neural networks. As will be appreciated, thefinal prediction is based upon a weighted average of the committee ofneural networks. In accordance with the present invention, analternative approach to developing a probabilistic model is based uponthe Bayesian evidence approximation, which allows development of aprobabilistic model without sampling and, thus, allows the use of asingle neural network model in the probabilistic formulation.Accordingly, in the present application, a probabilistic model isdeveloped using Bayesian techniques for the model shown in FIG. 5. Thisprobabilistic model may be developed without sampling. According tocertain embodiments, it is shown how this model may be trained. Inaddition, the present application shows how the quality of the model canbe computed using a confidence metric. This confidence metric may thenbe used to determine whether additional data is required for trainingthe model and what that additional data is. Also, it is shown how theprobabilistic model can practically be used in an optimization tocompute a set of MVs that can be used to control an industrial system,such as the plant of FIG. 1.

Training Algorithm for Probabilistic Neural Networks with DisturbanceRejection

In this section, a training algorithm for probabilistic neural networkswith disturbance rejection will be discuss that accords with exemplaryembodiments of the present invention. First, with reference to U.S. Pat.No. 7,123,971, a review of a training algorithm for deterministic neuralnetworks with disturbance rejection is provided. Following that, analgorithm for the calculation of the second derivative of an errorfunction with respect to the weights is provided. Then, the Bayesianapproach to the neural network model in FIG. 5 is derived, and anapproach to Bayesian hyperparameter re-estimation is shown that is inaccordance with the present invention. In the final section of this partof the application, a review is provided of an exemplary trainingalgorithm for the probabilistic neural network with disturbancerejection that accords with the present invention.

The disturbance rejection algorithm, which may be used to train theneural network model in FIG. 5, will now be discussed. It is describedin full detail in U.S. Pat. No. 7,123,971, which is incorporated hereinby reference in its entirety. First, the model to be trained is definedusing the disturbance rejection algorithm along with a training dataset.As shown in FIG. 3, the training dataset may be generated via collectingdata from time 0 to T for the correlated MVs, DVs and CVs of the system.The MVs and DVs at time (t) may be combined into a generic input vector,x_(t)=[x_(t,1) x_(t,2) . . . x_(t,X)]^(T), where X is the number ofinputs. For convenience and without loss of generality, a model with asingle output or CV is used for the purposes of this example. Thesampled data for the output of interest at time (t), which may be a CVsuch a NOx, is represented by the scalar d_(t). Thus, an exemplaryneural network model at time (t) and weights w may be defined as:

m _(t)=NN(x _(t−1) ,w)  Equation 1

where m=[m₁ m₂ . . . m . . . m_(T)]^(T) refers to the predicted outputvector from the model. The predicted output at time (t), y_(t), may bedefined as:

y _(t) =m _(t) +k(d _(t−1) −m _(t−1))  Equation 2

The vector y=[y₁ y₂ . . . y_(t) . . . y_(T)]^(T) is the predicted outputover the dataset. Using a neural network as the model results in thefollowing equation representing the model shown in FIG. 5:

y _(t)=NN(x _(t−1) ,w)+k(d _(t−1)−NN(x _(t−2) ,w))  Equation 3

Accordingly, as used herein, the following definitions are provided:

x_(t) is the input vector in the historical dataset, comprised of theMVs and DVs;

d_(t) is the output in the historical dataset;

d=[d₁ d₂ . . . d_(t) . . . d_(T)]^(T) is a vector contain all outputs inthe historical dataset;

m_(t) is the predicted output from the models;

m=[m₁ m₂ . . . m . . . m_(T)]^(T) is the predicted output vector fromthe models;

NN(x_(t),w) is a neural network prediction at time (t) given a set ofweights; and

w is the weight vector of a neural network.

To find the weight vector w of the neural network, the following errorfunction may be minimized:

$\begin{matrix}{E = {{\frac{1}{2}\alpha \; w^{T}w} + {\frac{1}{2}{\sum\limits_{t = 3}^{T}\; \left( {y_{t} - d_{t}} \right)^{2}}}}} & {{Equation}\mspace{14mu} 4}\end{matrix}$

where the hyperparameter a is a penalty on the size of the weights inthe weight vector w.

Alternatively, using:

ε_(t) =y _(t) −d _(t)  Equation 5

the error function can be written as:

$\begin{matrix}{E = {{\frac{1}{2}\alpha \; w^{T}w} + {\frac{1}{2}{\sum\limits_{t = 2}^{T}\varepsilon_{t}^{2}}}}} & {{Equation}\mspace{14mu} 6}\end{matrix}$

The Disturbance Rejection Algorithm

Given the definitions above, the disturbance rejection algorithm of U.S.Pat. No. 7,123,971 is provided below. As will be appreciated, thisalgorithm is based on a gradient descent approach or a modified versionof backpropagation. With reference again to U.S. Pat. No. 7,123,971, thedisturbance rejection algorithm may include the following steps and maybe employed to update the weights of the model at each trainingiteration:

-   -   i) Forward Iteration of the Model: A forward propagation of the        model across the time series, from (t=2) to (t=T) is performed        using the following equation:

y _(t)=NN(x _(t−1) ,w)+k(d _(t−1)−NN(x _(t−2) ,w))  Equation 7

-   -   ii) Modified Backpropagation: The derivative of the error with        respect to the weights and feedback coefficient are computed        using:

$\begin{matrix}{\frac{\partial E}{\partial w} = {{\alpha \; w} + {\sum\limits_{t = 1}^{T}{\left( {\varepsilon_{t} - {k\; \varepsilon_{t + 1}}} \right)\frac{\partial{{NN}\left( {x_{t - 1},w} \right)}}{\partial w}}}}} & {{Equation}\mspace{14mu} 8} \\{\frac{\partial E}{\partial k} = {\sum\limits_{t = 2}^{T}{{- \varepsilon_{t}}\varepsilon_{t - 1}^{\prime}}}} & {{Equation}\mspace{14mu} 9} \\{\varepsilon_{t} = {y_{t} - d_{t}}} & {{Equation}\mspace{14mu} 10} \\{\varepsilon_{t}^{\prime} = {{{NN}\left( {x_{t},w} \right)} - d_{t}}} & {{Equation}\mspace{14mu} 11}\end{matrix}$

-   -    Where ε_(T+1)=ε₁=0. The derivative of the error with respect to        the weights can be implemented by backpropagating ε_(t)−kε_(t+1)        through the network at each time interval from 2 to T−1. At time        T, ε_(t), is backpropagated through the network while at time 1,        −kε₂, is backpropagated.    -   iii) Model Parameter Update: The weights may then be updated        using steepest descent or other gradient based weight learning        algorithms. To guarantee stability of the time series model, it        generally is advisable to constrain the feedback coefficient, k,        to a value of less than 1.0 and greater than 0.0.

Matrix Version of the Disturbance Rejection Algorithm

For convenience in deriving the second derivative of the disturbancerejection model in FIG. 5, a matrix version of the disturbance rejectionalgorithm may be used, as provided below. To define the matrix of thedisturbance rejection algorithm, the following T×T matrices (T is thenumber of samples in the dataset) may be needed:

$\begin{matrix}{C = \begin{bmatrix}0 & 0 & 0 & \; & 0 & 0 & 0 \\0 & 1 & 0 & \ldots & 0 & 0 & 0 \\0 & 0 & 1 & \; & 0 & 0 & 0 \\\; & \vdots & \; & \ddots & \; & \vdots & \; \\0 & 0 & 0 & \; & 1 & 0 & 0 \\0 & 0 & 0 & \ldots & 0 & 1 & 0 \\0 & 0 & 0 & \; & 0 & 0 & 1\end{bmatrix}} & {{Equation}\mspace{14mu} 12} \\{K = \begin{bmatrix}0 & 0 & 0 & \; & 0 & 0 & 0 \\k & 0 & 0 & \ldots & 0 & 0 & 0 \\0 & k & 0 & \; & 0 & 0 & 0 \\\; & \vdots & \; & \ddots & \; & \vdots & \; \\0 & 0 & 0 & \; & 0 & 0 & 0 \\0 & 0 & 0 & \ldots & k & 0 & 0 \\0 & 0 & 0 & \; & 0 & k & 0\end{bmatrix}} & {{Equation}\mspace{14mu} 13} \\{S = \begin{bmatrix}0 & 0 & 0 & \; & 0 & 0 & 0 \\1 & 0 & 0 & \ldots & 0 & 0 & 0 \\0 & 1 & 0 & \; & 0 & 0 & 0 \\\; & \vdots & \; & \ddots & \; & \vdots & \; \\0 & 0 & 0 & \; & 0 & 0 & 0 \\0 & 0 & 0 & \ldots & 1 & 0 & 0 \\0 & 0 & 0 & \; & 0 & 1 & 0\end{bmatrix}} & {{Equation}\mspace{14mu} 14}\end{matrix}$

It should be noted that these matrices may be implemented using sparsematrices methods to allow efficient storage and calculation with them.Using these definitions of C and K, the forward propagation equation ofthe disturbance rejection algorithm, Equation 7, can be rewritten inmatrix form as:

y=C(m+K(d−m))  Equation 15

where:

m=[NN(x ₀ ,w)NN(x ₁ ,w), . . . ,NN(x _(T−1) ,w)]^(T)  Equation 16

Given ε=[ε₁ ε₂ . . . ε_(t) . . . ε_(T)]^(T), the error function can bewritten as:

E=½αw ^(T) w+½ε^(T) Cε  Equation 17

Then, the derivative of the error with respect to the weights, given byEquation 8, can be rewritten as:

$\begin{matrix}{\frac{\partial E}{\partial w} = {{\alpha \; w} + {\varepsilon^{T}{C\left( {I - K} \right)}\frac{\partial m}{\partial w}}}} & {{Equation}\mspace{14mu} 18}\end{matrix}$

where

$\frac{\partial m}{\partial w}$

is a T×N matrix, and T is the number of samples in the dataset, and N isthe number of weights, and I is the identity matrix of size T×T.

Similarly, the derivatives of the error with respect to the feedbackcoefficient of Equation 9 may be given by:

$\begin{matrix}{\frac{\partial E}{\partial k} = {{- \varepsilon^{T}}{CS}\; \varepsilon^{\prime}}} & {{Equation}\mspace{14mu} 19}\end{matrix}$

where ε′=[ε′₁ ε′₂ . . . ε′_(t) . . . ε′_(T)]^(T). Further, it can beshown that:

ε=C(I−K)ε′  Equation 20

Thus, the equations above may be rewritten as:

$\begin{matrix}{\frac{\partial E}{\partial w} = {{\alpha \; w} + {{{\varepsilon^{\prime}}^{T}\left( {I - K} \right)}^{T}{C\left( {I - K} \right)}\frac{\partial m}{\partial w}}}} & {{Equation}\mspace{14mu} 21} \\{and} & \; \\{\frac{\partial E}{\partial k} = {{- {{\varepsilon^{\prime}}^{T}\left( {I - K} \right)}^{T}}{CS}\; \varepsilon^{\prime}}} & {{Equation}\mspace{14mu} 22}\end{matrix}$

As will be appreciated, these equations are needed to compute theHessian, which is described in the next section.

Computing the Hessian

To implement the probabilistic version of the disturbance rejectionalgorithm, the data Hessian, which is the second derivative of the datacomponent of the error with respect to the weights, may be required. TheHessian, H, may be defined as:

$\begin{matrix}{H = \frac{\partial^{2}E_{D}}{\partial^{2}w}} & {{Equation}\mspace{14mu} 23}\end{matrix}$

where the data component of the error is defined by:

E _(D)=½ε^(T) Cε  Equation 24

Further, the Hessian may be computed using a modified version of theR-propagation algorithm. As will be appreciated, R-propagation is anefficient algorithm for multiplying a vector by the Hessian. The vector,v, may be defined as:

v=[v ₁ v ₂ . . . v _(t) . . . v _(N)]^(T)  Equation 25

where N is the number of weights in the weight vector. Multiplying v byH gives:

$\begin{matrix}{{v^{T}H} = {v^{T}\frac{\partial^{2}E_{D}}{\partial^{2}w}}} & {{Equation}\mspace{14mu} 26}\end{matrix}$

which can be written as:

$\begin{matrix}{{v^{T}H} = {v^{T}\frac{\partial}{\partial w}\left( \frac{\partial E_{D}}{\partial w} \right)}} & {{Equation}\mspace{14mu} 27}\end{matrix}$

which can be rewritten to:

$\begin{matrix}{{v^{T}H} = {R\left\{ \frac{\partial E_{D}}{\partial w} \right\}}} & {{Equation}\mspace{14mu} 28}\end{matrix}$

where R{ } is the R-propagation operator. The ith column of the Hessiancan be computed by selecting the vector v to be a vector of length Ncomposed of all zeros except the ith element which contains a 1. Byextension, the Hessian can be computed using N vector multiples toderive the N columns of the Hessian.

To compute the R-propagation value, the forward propagation needs to beexpressed through a neural network in matrix form. First, the weightmatrices used in the matrix form is defined. The N×1 weight vector ofthe neural network and its components are defined as:

w=[w ₁ b ₁ w ₂ b ₂]^(T)  Equation 29

where w₁ are the weights in the first layer and given by:

w ₁ =[w _(1,1) w _(1,2) . . . w _(1,H)]  Equation 30

where w_(1,1) are the weights between the inputs and the first hiddenunit and are given by:

w _(1,1) =[w _(1,1,1) w _(1,1,2) . . . w _(1,1,X)]  Equation 31

where w_(1,1,1) is a weight in the first layer between the first inputand first hidden unit. The weight vector between the input layer and theHth hidden unit is given by:

w _(1,H) =[w _(1,H,1) w _(1,H,2) . . . w _(1,H,X)]  Equation 32

where w_(1,H,N) is the weight in the first layer between the Nth inputand Hth hidden unit. (X is the number of inputs and H is the number ofinputs.)

The biases in the first layer may be defined as:

b ₁ =[b _(1,1) b _(1,2) . . . b _(1,H)]  Equation 33

where b_(1,1) is the bias to the first hidden unit. The second layerweight vector may defined as:

w ₂ =[w _(2,1) w _(2,2) . . . w _(2,H)]  Equation 34

where w_(2,1) is the weight between the first hidden unit and theoutput. Finally, b₂, is the bias of the output unit.

The input layer weight matrix may be defined as:

W ₁ =[w _(1,1) ^(T) w _(1,2) ^(T) . . . w _(1,H) ^(T)]  Equation 35

Finally, a flattening function which converts a matrix of weights in thefirst layer into a weight vector of the first layer may be defined as:

w ₁ =F(W ₁)  Equation 36

For convenience, a 1×N vector v is defined in a similar manner:

v=[v ₁ p ₁ v ₂ p ₂]^(T)  Equation 37

where v₁ has the same structure as the weights in the first layer andgiven by:

v ₁ =[v _(1,1) v _(1,2) . . . v _(1,H)]  Equation 38

where v_(1,1) is given by:

v _(1,1) =[v _(1,1,1) v _(1,1,2) . . . v _(1,1,X)]  Equation 39

and v_(1,H) is:

v _(1,H) =[v _(1,H,1) v _(1,H,2) . . . v _(1,H,X)]  Equation 40

The vector p₁ is given by:

p ₁ =[p _(1,1) p _(1,2) . . . p _(1,H)]  Equation 41

The v₂ vector is:

v ₂ =[v _(2,1) v _(2,2) . . . v _(2,H)]  Equation 42

Finally, p₂ is a scalar.

Furthermore, V₁ is defined as:

V ₁ =[v _(1,1) ^(T) v _(1,2) ^(T) . . . v _(1,H) ^(T)]  Equation 43

To write the matrix formulation, the inputs in vector form may bedefined as:

$\begin{matrix}{X = \begin{bmatrix}x_{0}^{T} \\x_{1}^{T} \\\vdots \\x_{T - 1}^{T}\end{bmatrix}} & {{Equation}\mspace{14mu} 44}\end{matrix}$

The matrix (or vector), O_(j,k), may be defined to be composed of amatrix with j rows and k columns and that contains ones in all elements.The hidden unit outputs over the dataset may be computed using:

A=XW ₁ +O _(T,1) b ₁  Equation 45

and

Z=h(A)  Equation 46

where h( ) is the hidden layer operator (typically the hyperbolictangent) that operates on each element of the matrix, A. Thus, Z,contains the output of the hidden units across all samples in thedataset. The output of the neural network over the dataset may be givenby:

m=Zw ₂ ^(T) +O _(T,1) b ₂  Equation 47

Finally, the last step of the forward propagation may be computing theoutput y which is given by Equation 15 and is shown again here forconvenience:

y=C(m+K(d−m))  Equation 48

Now the R-propagation values for A, Z, and m are shown. Starting with A:

R{A}=XV ₁ +O _(T,1) p ₁  Equation 49

To compute the R-propagation values of Z and other R-propagated values,the first and second derivatives of activation layer of the neuralnetwork may be needed. For the first derivative, the interest is in thefirst derivative of the output activation with respect to the input ofthe activation layer for each data sample. The first derivative of theoutput activations with respect to the input activations over thedataset can be computed by:

Z′=O _(T,N) +Z∘Z  Equation 50

where Z∘Z represents the Hadamard product (element multiplication).Similarly, for the second derivative, the interest in the secondderivative of the output activation with respect to the input of theactivation layer for each data sample. The second derivative of theoutput activations with respect to the input activations over thedataset can be computed by:

Z″=−2Z∘Z′  Equation 51

The R-propagation value of Z is given by:

R{Z}=Z′∘R{A}  Equation 52

and R{m} is given by:

R{m}=R{Z}w ₂ ^(T) +Zv ₁ +O _(T,1) p ₂  Equation 53

The value for R{y} can computed based on the definition of y given byEquation 48,

R{y}=C(I−K)R{m}  Equation 54

The next step is to define the backpropagation through the disturbancerejection layer into the neural network. Define u₂ to be the derivativeof error, E, with respect to the output of the neural network models, m.The vector u₂ is the backpropagated error to the output of the neuralmodels and can be found using Equation 18 to be:

u ₂=(I−K)^(T) Cε  Equation 55

As seen previously, a modified version of the output error, ε, isbackpropagated through the neural network. Now U₁ can be defined asbeing the derivative of error, E, with respect to the first layeractivations of the neural network models over the dataset, u₂. The valueof U₁ is:

U ₁ =Z′∘(u ₂ w ₂)  Equation 56

Given the backpropagation values, the R-propagation can be determinedfor u₂ and U₁. Using Equation 55, it is found:

R{u ₂}=(I−K)^(T) CR{y}  Equation 57

The value for R{U₁} is given as:

R{U ₁ }=Z″∘R{A}∘(u ₂ w ₂)+Z′∘(u ₂ v ₂)+Z′∘(R{u ₂ }w ₂)  Equation 58

The R-propagated values of the derivative of the error with respect tothe weights are given by:

$\begin{matrix}{{R\left\{ \frac{\partial E_{D}}{\partial w} \right\}} = \begin{bmatrix}{R\left\{ \frac{\partial E_{D}}{\partial w_{1}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial b_{1}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial w_{2}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial b_{2}} \right\}}\end{bmatrix}} & {{Equation}\mspace{14mu} 59}\end{matrix}$

where each of the terms of the equation are given by the following:

$\begin{matrix}{{R\left\{ \frac{\partial E_{D}}{\partial w_{1}} \right\}} = {F\left( {X^{T}R\left\{ U_{1} \right\}} \right)}} & {{Equation}\mspace{14mu} 60} \\{{R\left\{ \frac{\partial E_{D}}{\partial b_{1}} \right\}} = {O_{1,T}R\left\{ U_{1} \right\}}} & {{Equation}\mspace{14mu} 61} \\{{R\left\{ \frac{\partial E_{D}}{\partial w_{2}} \right\}} = {{Z^{T}R\left\{ u_{2} \right\}} + {R\left\{ Z \right\} u_{2}}}} & {{Equation}\mspace{14mu} 62} \\{{R\left\{ \frac{\partial E_{D}}{\partial b_{2}} \right\}} = {O_{1,T}R\left\{ u_{2} \right\}}} & {{Equation}\mspace{14mu} 63}\end{matrix}$

where F( ) is the flattening function used to convert from a matrix to avector. Thus, using Equation 28, it is found that:

$\begin{matrix}{{v^{T}H} = \begin{bmatrix}{R\left\{ \frac{\partial E_{D}}{\partial w_{1}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial b_{1}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial w_{2}} \right\}} & {R\left\{ \frac{\partial E_{D}}{\partial b_{2}} \right\}}\end{bmatrix}} & {{Equation}\mspace{14mu} 64}\end{matrix}$

where the terms of Equation 64 are defined by Equation 60-Equation 63.Once again, the full Hessian can be computed by doing N multiplicationof v^(T)H with the vector v^(T) varying in each multiplication such thatall elements are zero except for the ith element which is 1 in the ithmultiplication.

Thus, as just provided, it has been shown how to compute the datacomponent of the Hessian (the second derivative of the data component ofthe error with respect to the weights) given a weight vector w and adataset d and X. As just discussed, the calculation may includeEquations 44 through 63.

Disturbance Rejection with Additional Hyperparameters

To use Bayesian techniques on the disturbance rejection model of FIG. 5,it may be necessary to introduce more hyperparameters into thedefinition of the error function used in the disturbance rejectionalgorithm. The disturbance rejection algorithm uses only onehyperparameter in error, a. The matrix version of the error function wasgiven as:

E=½αw ^(T) w+½ε^(T) Cε  Equation 65

where the hyperparameter, α, is a penalty on the weights. Thehyperparameter, α, is typical found using cross-validation. A secondhyperparameter, β, may be introduced in the error equation as follows:

E=½αw ^(T) w+½βε^(T) Cε  Equation 66

Again, the hyperparameters, α and β, can be found using cross-validationtechniques.

Next, the weights may be penalized independently by using a vector, α,defined as:

α=[α₁ α₂ . . . α_(t) . . . α_(N)]  Equation 67

where N is the number of weights in the weight vector w. The errorfunction with individual weight penalties then may be defined as:

E=½w ^(T)(Diag(α))w+½βε^(T) Cε  Equation 68

where Diag(α) is an N×N matrix with a on the diagonal and zeros in allother elements. From this definition, it can be observed that the firstweight in vector w is penalized by α₁, the second weight by α₂, etc.Using Equation 18, the derivative of the error may be computed withrespect to the weights to be:

$\begin{matrix}{\frac{\partial E}{\partial w} = {{{{Diag}(\alpha)}w} + {\varepsilon^{T}{C\left( {I - K} \right)}\frac{\partial m}{\partial w}}}} & {{Equation}\mspace{14mu} 69}\end{matrix}$

Thus, instead of penalizing all weights independently, it may bepreferable to introduce structure into the penalty. It may be preferableto penalize all weights in the first layer of the neural network fromthe same input equally. In addition, it may be preferable to penalizeall biases in the first layer equally and all weights in the secondlayer together. The output bias may be individually penalized. Giventhese preferences, an a can be found with the following structure:

α=[α_(w1) α_(w1) . . . α_(w1) α_(b1) α_(w2) α_(b2)]  Equation 70

where α_(w1) is repeated H times (for the H hidden units) and is definedas:

α_(w1)=[α_(w1,1) α_(w1,2) . . . α_(w1,X)]  Equation 71

where α_(w1,1) is used to equally penalize the weight associated withthe first input to all hidden layer units equally. Similarly, α_(b1) isused to penalize all biases in the first layer equally by α_(b1) and isdefined by:

α_(b1)=α_(b1) O _(1,H)  Equation 72

And, similarly, α_(w2) is used to penalize all biases in the first layerequally by α_(w2) and is defined by:

α_(w2)=α_(w2) O _(1,H)  Equation 73

Finally, the scalar α_(b2) may be used to penalize the second layer biasof the neural network. Because the number of independent hyperparametersin a is equal to N+3, which is the number of inputs plus 3 and thenumber of inputs is typically greater than 10 in many real-worldapplications, cross validation typically is not practical for findingthe values of α and β. The next section provides a method for findingthe hyperparameters while also allowing us to calculate the variance ofthe output. Finally, the disturbance rejection algorithm in matrix formwith additional hyperparameters is derived in this section. Forconvenience, the weights and feedback coefficient that minimize theerror function found using this algorithm shall be denoted as w_(MAP)and k*.

Bayesian Prediction

It may be desired to predict the probabilistic distribution of thegeneralized output, τ_(t), given a model, set of input vectors, x_(t−1)and x_(t−2), and the output at time (t−1), d_(t−1). (It should be notedthat these general variables representing the output and inputs vectorat any arbitrary time rather than data specifically in a trainingdataset.) The Bayesian evidence framework may be used to derive thedesired probabilistic distribution, as described below.

First, it may be assumed that the conditional distribution of the outputwith respect to the inputs, p(τ_(t)|x_(t−1),x_(t−2)) is Gaussian with amean dependent upon the output of disturbance rejection based model andthe variance defined by the inverse variance of β:

p(τ_(t) |x _(t−1) ,x _(t−2) ,d _(t−1) ,w,β)=

(τ_(t) |y(x _(t−1) ,x _(t−2) ,d _(t−1) ,w),β⁻¹)  Equation 74

where

(τ_(t)|y(x_(t−1),x_(t−2),w),β⁻¹) represents the normal distribution withmean y(x_(t−1),x_(t−2),w) and variance β⁻¹. The normal distribution thenmay be defined as:

$\begin{matrix}{{\left( {\left. \tau_{t} \middle| {y\left( {x_{t - 1},x_{t - 2},d_{t - 1},w} \right)} \right.,\beta^{- 1}} \right)} = {\frac{1}{\sqrt{2{\pi\beta}^{- 1}}}e^{\frac{- {({\tau_{t} - {y{({x_{t - 1},x_{t - 2},d_{t - 1},w})}}})}^{2}}{2\beta^{- 1}}}}} & {{Equation}\mspace{14mu} 75}\end{matrix}$

It is also desirous to define the prior distribution of the weights tobe Gaussian and have the following form:

p(w|α)=

(w|O,Diag(α⁻¹))  Equation 76

Thus, the weights are zero mean with a variance defined by Diag(α⁻¹).

Assuming a set of T independent identically distributed samples ofd_(t), x_(t−1) and x_(t−2) defined by the dataset D={d₂ . . . d_(T)} andx₁ . . . x_(T), the likelihood function is given by:

$\begin{matrix}{{p\left( {\left. D \middle| w \right.,\beta} \right)} = {\prod\limits_{t = 2}^{T}{\left( {\left. \tau_{t} \middle| {y\left( {x_{t - 1},x_{t - 2},d_{t - 1},w} \right)} \right.,\beta^{- 1}} \right)}}} & {{Equation}\mspace{14mu} 77}\end{matrix}$

and the posterior distribution is given by:

p(w|D,α,β)=p(w|α)p(D|w,β)  Equation 78

The weights associated with the maximum of the posterior distribution,w_(MAP), can found by maximizing the log of the posterior which is givenby:

$\begin{matrix}{{\ln \; {p\left( {\left. w \middle| D \right.,\alpha,\beta} \right)}} = {{{- \frac{1}{2}}w\; {{Diag}(\alpha)}w^{T}} - {\frac{\beta}{2}{\sum\limits_{t = 2}^{T}\left( {y_{t} - d_{t}} \right)^{2}}}}} & {{Equation}\mspace{14mu} 79}\end{matrix}$

As will be appreciated, this is equivalent to:

ln p(w|D,α,β)=−½w ^(T)(Diag(α))w−½βε^(T) Cε  Equation 80

From the definition of the error function in Equation 68, it can befound that:

ln p(w|D,α,β)=−E  Equation 81

Thus, minimizing the error function of Equation 68 is equivalent tomaximizing the posterior distribution of Equation 78. Hence, thedisturbance rejection algorithm given above can be used to find the setof weights that maximize the posterior, w_(MAP). The resulting feedbackcoefficient from the disturbance rejection algorithm is denoted as k*.

Given the weights that maximize the posterior, w_(MAP), the marginalprobability of an output prediction may be determined by:

p(τ_(t) |D,x _(t−1) ,x _(t−2) ,d _(t−1) ,w,α,β)=

(τ_(t) |y(x _(t−1) ,x _(t−2) ,d _(t−1) ,w _(MAP)),σ²(x _(t−1) ,x _(t−)  Equation 82

Thus, the prediction is normally distributed aroundy(x_(t−1),x_(t−2),d_(t−1)w_(MAP)) with the variance defined as:

σ²(x _(t−1) ,x _(t−2) ,d _(t−1) ,w _(MAP))=β⁻¹ +g _(t) ^(T) A ⁻¹ g_(t)  Equation 83

where A is the second derivative of the output with respect to theweights and is given by:

A=Diag(α)I+βH  Equation 84

where the R-propagation algorithm output in the previous section can beused to compute the Hessian, H. The vector, g_(t) may be defined as:

$\begin{matrix}{{g_{t} = \frac{\partial{y\left( {x_{t - 1},x_{t - 2},w} \right)}}{\partial w}}}_{w = w_{MAP}} & {{Equation}\mspace{14mu} 85}\end{matrix}$

which using Equation 3 can be written as:

$\begin{matrix}{g_{t} = \left( {\frac{\partial{{NN}\left( {x_{t - 1},w} \right)}}{\partial w} - {k^{*}\frac{\partial{{NN}\left( {x_{t - 2},w} \right)}}{\partial w}}} \right._{w = w_{MAP}}} & {{Equation}\mspace{14mu} 86}\end{matrix}$

Since the gradient of the output of neural network with respect to theweights can be computed by backpropagating a value of 1 into thenetwork, the two partials of the equation above can be computed bybackpropagating a 1 into the networks given by NN(x_(t−1), w_(MAP)) andNN(x_(t−2),w_(MAP)). Using the definition of the normal distribution inEquation 82, it is found that the output distribution is:

$\begin{matrix}{{p\left( {\left. \tau_{t} \middle| D \right.,x_{t - 1},x_{t - 2},d_{t - 1},w,\alpha,\beta} \right)} = \frac{e^{\frac{- {({\tau_{t} - {y{({x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}})}}})}^{2}}{2{\sigma^{2}{({x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}})}}}}}{\sqrt{2{{\pi\sigma}^{2}\left( {x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}} \right)}}}} & {{Equation}\mspace{14mu} 87}\end{matrix}$

Given α, β and the training set, D, the disturbance rejection algorithmis used to find w_(MAP). Thus, given w_(MAP) along with thehyperparameters α, β, the probability distribution for the model of FIG.5 can be written as a function of only the inputs to the model:

$\begin{matrix}{{p\left( {\left. \tau_{t} \middle| x_{t - 1} \right.,x_{t - 2},d_{t - 1}} \right)} = \frac{e^{\frac{- {({\tau_{t} - {y{({x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}})}}})}^{2}}{2{\sigma^{2}{({x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}})}}}}}{\sqrt{2{{\pi\sigma}^{2}\left( {x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}} \right)}}}} & {{Equation}\mspace{14mu} 88}\end{matrix}$

Determining the Hyperparameters

By starting with an initial guess for α and β, which is typically 1 forall values, the weights that maximize the posterior distribution,w_(MAP), and the trained feedback coefficient, k*, can be determinedusing the disturbance rejection algorithm summarized above. Conventionalmethods provide algorithms for updating α and β given a dataset andw_(MAP) in the case of a standard neural network. In the presentapplication, these algorithms have been extended to the disturbancerejection based neural network model of FIG. 5.

Before presenting the updated algorithm, the vector q(α_(w1,1)) mustfirst be defined as:

q(α_(w1,1))=[q ₁(α_(w1,1)) . . . q _(i)(α_(w1,1)) . . . q_(N)(α_(w1,1))]  Equation 89

where α_(w1,1) is contained in the vector α as defined in Equation 70, Nis the number of weights, and q_(i)(α_(w1,1)) is defined as:

$\begin{matrix}{{q_{i}\left( \alpha_{{w\; 1},1} \right)} = \left\{ \begin{matrix}1 & {{{if}\mspace{14mu} \alpha_{i}} = \alpha_{{w\; 1},1}} \\0 & {{{if}\mspace{14mu} \alpha_{1}} \neq \alpha_{{w\; 1},1}}\end{matrix} \right.} & {{Equation}\mspace{14mu} 90}\end{matrix}$

where α_(i) is the ith element of the vector α. Thus, the vectorq(α_(w1,1)) denotes with a 1 where α_(w1,1) appears in the vector α.Similarly, the vectors q(α_(w1,i)), q(α_(w2)), and q(α_(b2)) can bedefined to denote where (α_(w1,i)), (α_(w2)), and (α_(b2)) appear in thevector α. For convenience, the set {q₁, . . . , q_(i), . . . , q_(X),q_(X+1), q_(X+2), q_(X+3)} is defined by the following:

q ₁ =q(α_(w1,1))  Equation 91

q _(i) =q(α_(w1,i))  Equation 92

q _(X) =q(α_(w1,X))  Equation 93

q _(X+1) =q(α_(b1))  Equation 94

q _(X+2) =q(α_(w2))  Equation 95

q _(X+3) =q(α_(b2))  Equation 96

where X is the number of inputs in the network. For convenience, thevector α is defined as:

α=[α_(w1,1) . . . α_(w1,i) . . . α_(w1,X) α_(b1) α_(w2)α_(w2)]  Equation 97

Furthermore, the scalar γ_(i) may be defined as:

$\begin{matrix}{\gamma_{i} = {{\sum\limits_{j = 1}^{N}q_{i,j}} - {a_{i}\left( {q_{i}\left( {{diag}\left( H^{- 1} \right)} \right)} \right.}}} & {{Equation}\mspace{14mu} 98}\end{matrix}$

where q_(i,j) is the jth element of the vector q_(i). Given thedefinition for the scalar γ_(i), the vector γ is:

γ=[γ₁ . . . γ_(i) . . . γ_(X+3)]  Equation 99

The updated versions of the elements of the vector α are given by:

$\begin{matrix}{a_{i} = \frac{\gamma_{i}}{\left( {q_{i} \circ w_{MAP}} \right)w_{MAP}^{T}}} & {{Equation}\mspace{14mu} 100}\end{matrix}$

An updated version of all X+3 elements of the vector α can be computedand the resulting values can be used to update the values of the vectorα.

The update equation then for β is given by:

$\begin{matrix}{\beta = \frac{N - {\sum\limits_{i = 1}^{X + 3}\gamma_{i}}}{\varepsilon^{T}C\; \varepsilon}} & {{Equation}\mspace{14mu} 101}\end{matrix}$

Then, using Equation 100 and Equation 101, the hyperparameters α and βcan be updated. It should be noted that the re-estimation given inEquation 100 and Equation 101 are often repeated N_(b) times in a row tore-estimated α and β where N_(b) is a user specified number, which was 3in the prior exemplary implementation.

Bayesian Disturbance Rejection Algorithm

FIG. 6 shows a flow diagram of the Bayesian Disturbance RejectionAlgorithm which is used to compute the trained values of thehyperparameters α and β, the weights, w_(MAP), and the feedbackcoefficient, k* of the model in FIG. 5. The trained values cansubsequently be used in to compute the probability distribution of theoutput of FIG. 5 using the Bayesian Disturbance Rejection PredictionCalculation.

To use the algorithm of FIG. 6, the following is needed: a trainingdataset d and X (typically collected using Design of Experiment), aninitial set of weights w (typically randomly generated), an initialvalue of k (typically 1), an initial set of hyperparameters α and β(typically all values are set to 1), a user specified number of times torepeat the main loop, N_(D)(typically 5), and a user specified number oftimes to repeat the re-estimation, N_(b)(typically 3).

Given the initial values, the Disturbance Rejection Algorithm isexecuted first as shown in FIG. 6. The results are then used to computethe data component of the Hessian using the Hessian calculation outlinedabove. Given the results of the Disturbance Rejection Algorithm and theHessian Calculation, the hyperparameters α and β are re-estimated usingthe Hyperparameter Re-estimation calculation. Next, a counter of thenumber of time through the main loop, i, is updated (i is initialized at0). If the main loop counter, i, is less than N_(D) the main loop isrepeated otherwise the algorithm exits with the trained values for thehyperparameters α and β, the weights, w_(MAP), and the feedbackcoefficient, k* of the model.

The algorithm shown in FIG. 6 is the preferred embodiment. However,extensions may be used which rely on cross-validation techniques for theselection of the following initial values: hyperparameters α and β,number of times to repeat the main loop, N_(D)(typically 5), and numberof times to repeat the re-estimation, N_(b)(typically 3). Specifically,the hyperparameters α and β have be initialized based upon the resultsof a cross-validation of initial hyperparameters over the set of values{0.001, 0.01, 0.1, 1, 10, 100}. Using cross-validation to determine theinitial hyperparameters for the algorithm of FIG. 6 has proved to beuseful.

Calculating Overall Model Confidence

Once the model is trained, the model of the form shown in FIG. 5 is usedfor optimization. Often, the mean value of the prediction is used whichis given by y(x_(t−1),x_(t−2),d_(t−1),w_(MAP)). To perform optimizationcorrectly, it is desirable to be able to predict the correct sign of thedirection of the change in the output due to changes in the input.Hence, it is desirable to know for a trained model what the probabilitythat the sign of the predicted change in the output is correct. Thisprobability can be defined as:

$\begin{matrix}{p\left( {\frac{{sign}\left( {\tau_{t} - d_{t - 1}} \right)}{{sign}\left( {{y\left( {x_{t - 1},x_{t - 2},d_{t - 1},w_{MAP}} \right)} - d_{t - 1}} \right)} > 0} \right)} & {{Equation}\mspace{14mu} 102}\end{matrix}$

where the function sign( ) is used to compute the sign of the values.

The probability distribution above can be approximated by sampling theprobability distribution over the training dataset. The probabilitysample for a given point in the dataset may computed as shown in FIG. 7.The probability distribution computed for any given point in the datasetcan be computed using the Bayesian Disturbance Rejection PredictionCalculation. The resulting probability distribution at time (t) is shownas the normal distribution shown in FIG. 7 with a mean ofy(x_(t−1),x_(t−2),d_(t−1),w_(MAP)) and a variance ofσ²(x_(t−1),x_(t−2),d_(t−1),w_(MAP))). In the example in FIG. 7, thevalue of sign(y(x_(t−1),x_(t−2),d_(t−1),w_(MAP))−d_(t−1)) is positive.The probability that sign(τ_(t)−d_(t−1)) is also positive is given bythe area under the normal distribution that is greater than d_(t−1).

The area under the normal distribution can be computed using the errorfunction. The integration of a Gaussian center at zero with a standarddeviation of σ is given by:

$\begin{matrix}{{\Phi (x)} = {\frac{1}{2} + {\frac{1}{2}{{erf}\left( \frac{x}{\sqrt{2}\sigma} \right)}}}} & {{Equation}\mspace{14mu} 103}\end{matrix}$

where erf( ) is the Gaussian error function. Given the known output attime t−1 which is given by d_(t−1), the probability that the predictionof the output at time t which is given by Equation 82 has the correctsign is given by:

$\begin{matrix}{{\Phi \left( {{{y\left( {x_{t - 1},x_{t - 2},w_{MAP}} \right)} - d_{t - 1}}} \right)} = {\frac{1}{2} + {\frac{1}{2}{{erf}\left( \frac{{{y\left( {x_{t - 1},x_{t - 2},w_{MAP}} \right)} - d_{t - 1}}}{\sqrt{2}{\sigma^{2}\left( {x_{t - 1},x_{t - 2},w_{MAP}} \right)}} \right)}}}} & {{Equation}\mspace{14mu} 104}\end{matrix}$

The equation above gives the sampled probability at one data point inthe dataset. The overall sampled probability can be calculated byaveraging over the dataset. Furthermore, for convenience, reference willbe made to probability multiplying by 100% as the over model confidence.Thus, the overall confidence of the model, η, is defined as the mean ofthe confidence over a dataset (typically, the training set) and is givenby:

$\begin{matrix}{\eta = {100\frac{\sum\limits_{t = 2}^{T}{C_{t,t}{\Phi \left( {{{y\left( {x_{t - 1},x_{t - 2},w_{MAP}} \right)} - d_{t - 1}}} \right)}}}{\sum\limits_{t = 2}^{T}C_{t,t}}}} & {{Equation}\mspace{14mu} 105}\end{matrix}$

where C_(t,t) is the t^(th) element of the diagonal of the matrix C.

Calculating Model Confidence Associated with Each Input

As will be appreciated, interest is not limited to overall confidence,but may also extends to the confidence associated with each of theinputs to the model, which is discussed in the following section. Tocomputed the confidence associated with each input, the input matrix, X,is expanded such that it contains individual input moves per variable.First, all moves in the input matrix may be identified by defining thematrix, B, which is a T−1×T matrix given by:

$\begin{matrix}{B = \begin{bmatrix}1 & {- 1} & 0 & \; & 0 & 0 & 0 \\0 & 1 & {- 1} & \ldots & 0 & 0 & 0 \\0 & 0 & 1 & \; & 0 & 0 & 0 \\\; & \vdots & \; & \ddots & \; & \vdots & \; \\0 & 0 & 0 & \; & {- 1} & 0 & 0 \\0 & 0 & 0 & \ldots & 1 & {- 1} & 0 \\0 & 0 & 0 & \; & 0 & 1 & {- 1}\end{bmatrix}} & {{Equation}\mspace{14mu} 106}\end{matrix}$

Using matrix B, Ω may be defined, which is a T−1×X matrix of thedifference in the inputs given by:

Ω=BX  Equation 107

Thus, Ω, contains the changes in the inputs over the dataset.Furthermore, the matrix, δ, may be defined which is a T×X matrix withelements given by:

$\begin{matrix}{\delta_{i,j} = \left\{ \begin{matrix}1 & {{{if}\mspace{14mu} \Omega_{i,j}} \neq 0} \\0 & {{{if}\mspace{14mu} \Omega_{i,j}} = 0}\end{matrix} \right.} & {{Equation}\mspace{14mu} 108}\end{matrix}$

Thus, the matrix, δ, indicates where moves in the inputs have occurredin the dataset.

To compute the confidence per inputs, a new input matrix, Γ, needs to bedefined as:

$\begin{matrix}{\Gamma = \begin{bmatrix}x_{0}^{T} \\x_{0}^{T} \\x_{1}^{T} \\x_{1}^{T} \\\vdots \\x_{T - 2}^{T} \\x_{T - 2}^{T}\end{bmatrix}} & {{Equation}\mspace{14mu} 109}\end{matrix}$

Furthermore, the vector, l_(j), needs to be defined as:

$\begin{matrix}{l_{j} = \begin{bmatrix}x_{0,j}^{T} \\x_{1,j}^{T} \\x_{1,j}^{T} \\x_{2,j}^{T} \\\vdots \\x_{{T - 2},j}^{T} \\x_{{T - 1},j}^{T}\end{bmatrix}} & {{Equation}\mspace{14mu} 110}\end{matrix}$

Finally, matrix, Λ^(j), needs to be defined as:

Λ^(j)=[Γ₁ Γ₂ . . . l _(j) . . . Γ_(N)]  Equation 111

where Γ_(i) represents the ith column of matrix Γ_(i). As will beappreciated, the matrix Λ^(j) contains the T-paired moves of the jthinput while all other variables are held equal. This matrix can then beused to compute the confidence for the jth input, ψ_(j), using

$\begin{matrix}{\psi_{j} = {100\frac{\sum\limits_{t = 2}^{T}{\delta_{{t - 1},j}C_{t,t}{\Phi \left( {{{y\left( {\Lambda_{2{({t - 1})}}^{j},\Lambda_{{2{({t - 1})}} - 1}^{j},w_{MAP}} \right)} - d_{t - 1}}} \right)}}}{\sum\limits_{t = 2}^{T}{\delta_{{t - 1},j}C_{t,t}}}}} & {{Equation}\mspace{14mu} 112}\end{matrix}$

where Λ^(j) _(i) is the ith row of the matrix Λ^(j). The confidence ofeach input is simply the average confidence of every move associate withthe input in the dataset while holding all other inputs constant. Theconfidence vector of all N inputs can be defined as:

ψ=[ψ₁ . . . ψ_(j) . . . ψ_(N)]  Equation 113

Design of Experiments Based on Model Confidence

As discussed in relation to FIG. 3, a “Design of Experiment” process maybe used to collect data for training the disturbance rejection model ofFIG. 5. FIG. 8 shows how the model confidence by input can be used toautomatically expand the collection of data used to train a model. Theflow diagram in FIG. 8 starts with the initial collection of data asshown in FIG. 3. In this case, the training data d and X which contain Ttraining exemplars is shown. Given the initial training data, as shownin FIG. 8, the Bayesian Disturbance Rejection algorithm can be used totrain the model of FIG. 5. Next, in FIG. 8, the Model Confidence byInput Calculation is performed. Results of the Confidence Calculationcan be compared to a threshold, in this case 90%, as shown in the tableillustrated in FIG. 9 for six MVs, MV₁-MV₆. If the confidence is abovethe threshold, no further additional data needs to be collected.However, if the data is below the acceptable confidence, additionalDesign of Experiments can be run for these MVs. For example, 10additional experiments could be performed for each MV that does notreach the acceptable confidence. As shown in FIG. 8, the additional datacan be added to the original training dataset to create a new dataset.Using this dataset, the model can be retrained. This process can berepeated until the necessary confidence per input is achieved. Thus, anautomated method for training a model can be implemented based upon theconfidence calculation.

An extension to the method shown in FIG. 8 allows automatic selection ofthe inputs based upon the confidence. It may be found that collectingadditional data in FIG. 8 does not results in the confidence exceedingthe threshold after collecting additional data. In this case, instead ofcontinuing to run the main loop of FIG. 8, it may be preferable to runthe main loop at most 3 times and then eliminate any inputs from themodel that have not achieved the desired confidence threshold. Once theinputs have been eliminated, the Bayesian Disturbance RejectionAlgorithm can be run to train the model and a final check can be done toverify that the remaining inputs are above the threshold. If not abovethe threshold, the inputs can be removed, training done again, and theconfidence threshold checked. This cycle can continue until all inputshave a confidence above the acceptable confidence.

The flow diagram of FIG. 10 may be used to do a daily update of themodel where additional data is collected each day for use in trainingthe model. Again, an initial set of data is used to train a model usingthe Bayesian Disturbance Rejection Algorithm of FIG. 6. Next, the modelconfidence per input is computed. At this point, a delay is introducedinto the loop typically of approximately one day. During this time, asdiscussed below, the model is used in the closed loop optimizer of FIG.2. However, after a 24 hour wait period, the delay is stopped andadditional data is collected as shown in FIG. 10. It should be notedthat the delay may be extended if other conditions are not met such asthe unit being at steady state. The confidence calculation per input isused in collecting additional data. Typically, two inputs are selectedrandomly based upon their confidence for design of experiments testingover a one hour period. If the sample period is 15 minutes, two movesper selected input are made during the data collection. The criteria forselecting which inputs to test is based upon the random draw of twoinputs where the probability of selection is a weighted inverse of theconfidence of the input. Thus, inputs with lower confidence are morelikely to be tested. As shown in FIG. 10, the additional collected datais added to the original data and the model is retrained. The main loopexecutes approximately once a day and data is added to the model fortraining based upon the model confidence by input calculation. It shouldbe noted that the Bayesian Disturbance Rejection Algorithm shown in FIG.10 could be replaced by the algorithm shown in FIG. 8. In addition, theinput selection algorithm described above could be used to the replacethe Bayesian Disturbance Rejection.

Optimizer

As shown in FIG. 2, an optimizer is used to minimize a “cost function”subject to a set of constraints. The cost function is a mathematicalrepresentation of a desired goal or goals. For instance, to minimizeNOx, the cost function includes a term that decreases as the level ofNOx decreases. One common method for minimizing a cost function is knownas “gradient descent optimization.” Gradient descent is an optimizationalgorithm that approaches a local minimum of a function by taking stepsproportional to the negative of the gradient (or the approximategradient) of the function at the current point.

Constraints may be placed upon both the inputs (MVs) and outputs (CVs)of the boiler at a future time. Typically, constraints that areconsistent with limits associated with the DCS are placed upon themanipulated variables. Constraints on the outputs (CVs) are determinedby the problem that is being solved.

A nonlinear model can be used to determine the relationship between theinputs and outputs of a boiler. Accordingly, a nonlinear programmingoptimizer is used to solve the optimization problem in accordance withthis embodiment of the present invention. However, it should beunderstood that a number of different optimization techniques may beused depending on the form of the model and the costs and constraints.For example, it is contemplated that the present invention may beimplemented by using, individually or in combination, a variety ofdifferent types of optimization approaches. These optimizationapproaches include, but not limited to, linear programming, quadraticprogramming, mixed integer non-linear programming (NLP), stochasticprogramming, global non-linear programming, genetic algorithms, andparticle/swarm techniques. See R. Baldick, Applied Optimization:Formulation and Algorithms for Engineering Systems, Cambridge UniversityPress, Cambridge, UK, 2009.

Given the cost function and constraints, a non-linear program (NLP)optimizer typically solves problems with 20 manipulated variables and 10controlled variables in less than one second. This is sufficiently fastfor most applications since the optimization cycle is typically in therange of anywhere from 20 seconds to 20 minutes. More details on theformulation of the cost function and constraints are provided in thereference S. Piche, B. Sayyar-Rodsari, D. Johnson and M. Gerules,“Nonlinear model predictive control using neural networks,” IEEE ControlSystems Magazine, vol. 20, no. 2, pp. 53-62, 2000, which is fullyincorporated herein by reference.

With reference now to FIG. 11, a form of the disturbance rejection modelis shown that is useful for optimization. In this case, it is desirableto determine the MVs at time (t) in order to minimize a cost functionwhich is based upon the CVs at time (t+1). In addition, the MVs and DVsat time (t−1) and the resulting CV at time (t) are known, which may beused to calculate a bias adjusted value k(NO_(x)(t)−NO_(x) ^(NN)(t))that is added to the neural network prediction at time (t+1), NO_(x)^(NN)(t+1). In the optimizer, the bias adjusted value,k(NO_(x)(t)−NO_(x) ^(NN)(t)), is fixed. In addition, the measureddisturbance value, DV(t), is fixed. Thus, only the values of MV(t) maybe adjusted by the optimizer to determine the CV of interested, in thiscase, NO_(x) ^(P)(t+1). Typically, the optimizer is executed at the samefrequency as the sampling period of the original data. If the samplingperiod was 15 minutes as shown in FIG. 3, then the optimizer would beexecuted every 15 minutes and the results of the optimizer would beoutput to the DCS and held for 15 minutes.

Optimization Example Using Mean Value of the Bayesian DisturbanceRejection Prediction

In this section, it is shown how a model of the form shown in FIG. 11which is trained using the Bayesian Disturbance Rejection Algorithm canbe used in an optimization solution to determine a set of MVs at time(t). Specifically, the mean value of the Bayesian Disturbance RejectionPrediction Calculation will be used in the optimization problem in thissection. In the next section, the probability distribution of theprediction is used in the optimization problem.

Without a loss of generality, the vector of MVs and DVs at time (t) canbe represented by the vector x_(t). The optimizer will be used todetermine the optimal values of the MVs, hence, only the terms of vectorx_(t) representing the MVs will be optimized while the other termsrepresenting the DVs will be held constant.

As shown in FIG. 11, in optimization, it is desirable to determine theoptimal values for the vector x_(t) at current time t to affect theoutput in the future y_(t+1) given the inputs from a previous time,x_(t−1), and the reading of the actual value of the output at thecurrent time, d_(t). Using Equation 2, it is found that:

y _(t+1) =m _(t+1) +k(d _(t) −m _(t))  Equation 114

Substituting in a neural network with trained weights, w_(MAP), andfeedback coefficient, k*, computed using the Bayesian DisturbanceRejection Algorithm (FIG. 6) or the algorithms shown in FIGS. 8 and 9,the following prediction of the output is provided:

y _(t+1)=NN(x _(t) ,w _(MAP))+k*(d _(t)−NN(x _(t−1) ,w_(MAP)))  Equation 115

It may be desirous to predict the value of NOx at time (t+1). Using thetraining algorithms described above, an optimal set of weights topredict NOx, which are w_(MAP,NOx) and the feedback coefficient,k*_(NOx), can be determine. The prediction of mean value of the NOx isgiven by:

y _(t+1,NOx)=NN(x _(t) ,w _(MAP,NOx))+k* _(NOx)(d _(t)−NN(x _(t−1) ,w_(MAP,NOx)))  Equation 116

Similarly, a prediction of CO could be defined as:

y _(t+1,CO)=NN(x _(t) ,w _(MAP,CO))+k* _(CO)(d _(t)−NN(x _(t−1) ,w_(MAP,CO)))  Equation 117

As shown above, optimization may be performed by defining a costfunction and a set of constraints and then minimizing the cost functionwhile observing the constraints. The cost function and constraintsrepresent the problem that one desires to solve. In an exemplary case,it is desirable to minimize the NOx, maintain the CO below a certainlevel, minimize the change in the inputs (MVs), and maintain the inputsin a bounded region. The cost function is used to express the desire tominimize the NOx, maintain the CO below a certain level, and minimizethe change in the inputs. Constraints are used to maintain the inputs ina bounded region.

An exemplary cost function to achieve these goals may be defined as:

$\begin{matrix}{{J\left( x_{t} \right)} = {{\frac{\pi_{1}}{2}y_{{t + 1},{NOx}}^{2}} + {\pi_{2}s\; e^{\frac{({y_{{t + 1},{CO}} - b_{co}})}{s}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)^{T}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 118}\end{matrix}$

where π₁, π₂, and π₃ are tuning constants for the optimization problemand are set by the user to achieve an appropriate desired balancebetween the three goals in the cost function. The first term of theequation above is used to minimize NOx, the second term is used topenalize CO over the bound b_(co) (s is used to scale the penalty), andthe final term is used to penalize movement in the inputs. Theconstraints are implement by the following bounds on the inputs:

b _(x,−) ≦x _(t) ≦b _(x,+)  Equation 119

where the vectors b_(x,−) and b_(x,+) represent the lower and upperbounds respectively on the input at time t.

According to exemplary embodiments, a gradient descent optimizer may beused to determine the minimal value of J(x_(t)) while observing theconstraints on the inputs. To perform gradient descent, the gradient ofthe cost function with respect to the inputs at time (t) is needed:

$\begin{matrix}{\frac{\partial{J\left( x_{t} \right)}}{\partial x_{t}} = {{\pi_{1}y_{{t + 1},{NOx}}\frac{\partial y_{{t + 1},{NOx}}}{\partial x_{t}}} + {\pi_{2}e^{\frac{({y_{{t + 1},{CO}} - b_{co}})}{s}}\frac{\partial y_{{t + 1},{CO}}}{\partial x_{t}}} + {\pi_{3}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 120}\end{matrix}$

Using Equation 116 and Equation 117, the previous equation can bewritten as:

$\begin{matrix}{\frac{\partial{J\left( x_{t} \right)}}{\partial x_{t}} = {{\pi_{1}y_{{t + 1},{NOx}}\frac{\partial{{NN}\left( {x_{t},w_{{MAP},{NOx}}} \right)}}{\partial x_{t}}} + {\pi_{2}e^{\frac{({y_{{t + 1},{CO}} - b_{co}})}{s}}\frac{\partial{{NN}\left( {x_{t},w_{{MAP},{CO}}} \right)}}{\partial x_{t}}} + {\pi_{3}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 121}\end{matrix}$

As will be shown, the first term of the previous equation can becomputed by backpropagating the value π₁y_(t+1,NOx) through the neuralnetwork represented NN(x_(t),w_(MAP,NOx)) to the inputs. Similarly, thesecond term can be computed by backpropagating the value

$\pi_{2}e^{\frac{({y_{{t + 1},{CO}} - b_{co}})}{s}}$

through the neural network represented NN(x_(t),w_(MAP,CO)) to theinputs. Thus, the derivative can be computed and steepest descentoptimization may be performed to calculated the optimal value of x_(t).The resulting optimized value for x_(t) is output to the DCS and heldconstant for 15 minutes at which point another optimization cycle isexecuted.

Optimization Using the Probability Distribution of the BayesianDisturbance Rejection Prediction

As will be appreciated, the previous section uses the mean value of theBayesian disturbance rejection prediction calculation in the costfunction of the optimization problem. In this section, the probabilisticform given by Equation 87 will be used in the cost function. The modelsused in this section were trained using the Bayesian disturbancerejection algorithm, of FIG. 6, or the algorithms shown in FIGS. 8 and9.

The probability distribution of the NOx at time t+1 is defined asp(τ_(t+1,NOx)|D_(NOx),x_(t),x_(t−1),w,β_(NOx)) is found using BayesianDisturbance Rejection algorithm given a training dataset D_(NOx). Thedistribution is given by:

p(τ_(t+1,NOx) |D _(NOx) ,x _(t) ,x _(t−1) ,w,β _(NOx))=

(τ_(t+1) |y(x _(t) ,x _(t−1) ,d _(t,NOx) ,k* _(NOx) ,w_(MAP,NOx)),σ_(NOx) ²(x _(t) ,x _(t−1) ,d _(t−1,NOx) ,k* _(NOx) ,w_(MAP,NOx).  Equation 122

where y(x_(t),x_(t−1),d_(t,NOx),w_(MAP,NOx)) is given by Equation 116and σ_(NOx) ²(x_(t),x_(t−1),d_(t,NOx),w_(MAP,NOx)) is given by:

σ_(NOx) ²(x _(t) ,x _(t−1) ,d _(t,NOx) ,k* _(NOx) ,w _(MAP))=β_(NOx) ⁻¹+g _(t+1,NOx) ^(T) A _(NOx) ⁻¹ g _(t+1,NOx)  Equation 123

where A_(NOx) is the second derivative of the output with respect to theweights and is given by:

A _(NOx)=Diag(α_(NOx))I+β _(NOx) H _(NOx)  Equation 124

where the R-propagation algorithm given above can be used to compute theHessian, H. The vector, g_(t+1), is defined as:

$\begin{matrix}{{g_{{t + 1},{NOx}} = \frac{\partial{y\left( {x_{t},x_{t - 1},d_{t,{NOx}},k_{NOx}^{*},w} \right)}}{\partial w}}}_{w = w_{{MAP},{NOx}}} & {{Equation}\mspace{14mu} 125}\end{matrix}$

Since the only set of variables that will be optimized using theoptimizer are the inputs at time (t), the mean and variance of Equation123 may be written as:

y _(MAP,NOx)(x _(t))=y(x _(t) ,x _(t−1) ,d _(t,NOx) ,k* _(NOx) ,w_(MAP,NOx))  Equation 126

and:

σ_(MAP,NOx) ²(x _(t))=σ_(NOx) ²(x _(t) ,x _(t−1) ,d _(t,NOx) ,k* _(NOx),w _(MAP))  Equation 127

Thus, the distribution for NOx is conveniently given by:

p(τ_(t+1,NOx) |D _(NOx) ,x _(t) ,x _(t−1) ,w,β _(NOx))=

(τ_(t+1) |y _(MAP,NOx)(x _(t)),σ_(MAP,NOx) ²(x _(t)))  Equation 128

Similarly, the distribution for CO is given by:

p(τ_(t+1,CO) |D _(CO) ,x _(t) ,x _(t−1) ,w,β _(CO))=

(τ_(t+1) |y _(MAP,CO)(x _(t)),σ_(MAP,CO) ²(x _(t)))  Equation 129

Once again, it may be desirous to minimize the NOx, maintain the CObelow a certain level, minimize the change in the inputs and maintainthe inputs in a bounded region. The cost function is used to express thedesire to minimize the NOx, maintain the CO below a certain level, andminimize the change in the inputs. The constraints are used to maintainthe inputs in a bounded region. Again, the optimization problem isdefined by a cost function and constraints. The cost function can bedefined as the following:

$\begin{matrix}{{J\left( x_{t} \right)} = {{\frac{\pi_{1}}{2}{\int_{- \infty}^{\infty}{\tau_{{t + 1},{NOx}}^{2}{p\left( {\left. \tau_{{t + 1},{NOx}} \middle| D_{NOx} \right.,x_{t},x_{t - 1},w,\beta_{NOx}} \right)}d\; \tau_{{t + 1},{NOx}}}}} + {\pi_{2}{\int_{- \infty}^{\infty}{s\; e^{\frac{({\tau_{{t + 1},{CO}} - b_{co}})}{s}}{p\left( {\left. \tau_{{t + 1},{NOx}} \middle| D_{NOx} \right.,x_{t},x_{t - 1},w,\beta_{NOx}} \right)}d\; \tau_{{t + 1},{CO}}}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)^{T}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 130}\end{matrix}$

The constraints are implement by the following bounds on the inputs:

b _(x,−) ≦x _(t) ≦b _(x,+)  Equation 131

As will be appreciated, the cost function and bound constraints are thesame as those defined in Equation 118 and Equation 119 except for thefirst two terms of Equation 130 are the integral of the penalty over thedistribution of the output. Substituting the normal distribution ofEquation 128 and Equation 129 into the cost function of Equation 130,the following is shown:

$\begin{matrix}{{J\left( x_{t} \right)} = {{\frac{\pi_{1}}{2}{\int_{- \infty}^{\infty}{\tau_{{t + 1},{NOx}}^{2}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}e^{- \frac{{({\tau_{t} - {y_{{MAP},{NOx}}{(x_{t})}}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}d\; \tau_{{t + 1},{NOx}}}}} + {\pi_{2}{\int_{- \infty}^{\infty}{s\; e^{\frac{({\tau_{{t + 1},{CO}} - b_{co}})}{s}}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}e^{- \frac{{({\tau_{t} - {y_{{MAP},{CO}}{(x_{t})}}})}^{2}}{2{\sigma_{{MAP},{CO}}^{2}{(x_{t})}}}}d\; \tau_{{t + 1},{CO}}}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)^{T}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 132}\end{matrix}$

An analytic solution to the cost function and the derivative of thefunction are not known. For this reason, an approximation to the costfunction is used:

$\begin{matrix}{{J\left( x_{t} \right)} \approx {{\frac{\pi_{1}}{2}{\sum\limits_{i = {- \infty}}^{\infty}{\left( {i\; \Delta} \right)^{2}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({{i\; \Delta} - {y_{{MAP},{NOx}}{(x_{t})}}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}\Delta}}} + {\pi_{2}{\sum\limits_{i = {- \infty}}^{\infty}{s\; e^{\frac{({{i\; \Delta} - b_{co}})}{s}}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({{i\; \Delta} - {y_{{MAP},{CO}}{(x_{t})}}})}^{2}}{2{\sigma_{{MAP},{CO}}^{2}{(x_{t})}}}}\Delta}}} + \; {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)^{T}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 133}\end{matrix}$

where the sum is equal to the integral as Δ goes to zero. To practicallycompute the summation, appropriate bounds for the summations and for Δin the equation above must be selected. Because the normal distributionis essentially 0 three standard deviation from the mean, it is requiredonly to sum around the mean value of the distribution within threestandard deviations. Furthermore, the values for Δ in each of sums inthe equation above can be defined based upon the standard deviations ofthe distribution. For this reason, Δ_(NOx) may be selected to be:

$\begin{matrix}{\Delta_{NOx} = \frac{6{\sigma_{{MAP},{NOx}}\left( x_{t} \right)}}{1000}} & {{Equation}\mspace{14mu} 134}\end{matrix}$

Similar Δ_(CO) can be defined for the CO variable. Given Δ_(NOx) andΔ_(CO), Equation 133 may be reasonably approximate by:

$\begin{matrix}{{J\left( x_{t} \right)} \approx {{\frac{\pi_{1}}{2}{\sum\limits_{i = {- 500}}^{500}{\left( {{i\; \Delta_{NOx}} + {y_{{MAP},{NOx}}\left( x_{t} \right)}} \right)^{2}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({i\; \Delta_{NOx}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}\Delta_{NOx}}}} + {\frac{\pi_{2}}{2}{\sum\limits_{i = {- 500}}^{500}{s\; e^{\frac{({{i\; \Delta_{CO}} + {y_{{MAP},{CO}}{(x_{t})}} - b_{co}})}{s}}\frac{1}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({i\; \Delta_{CO}})}^{2}}{2\; {\sigma_{{MAP},{CO}}^{2}{(x_{t})}}}}\Delta_{CO}}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)^{T}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 135}\end{matrix}$

As will be appreciated, the equation above provides a reasonable andpractical approach to computing the cost function in real time. To solvefor the gradient of the cost function with respect to the inputs at time(t), it is necessary to use an approximation. The easiest way toapproximate the gradient is to use numerical differential. However, thisapproach is not computation efficient and a better approach is sought.

Alternatively, it is observed that the mean values used in the equationabove for the cost function are significantly more affected by changesin the inputs than the variances. Thus, it is a reasonable approximationto assume the variances in the equation above are not affected bychanges in the input. Assuming constant variances for σ_(MAP,NOx)²(x_(t)) and σ_(MAP,CO) ²(x_(t)), the gradient of the cost function withrespect to the inputs is approximated by:

$\begin{matrix}{\frac{\partial{J\left( x_{t} \right)}}{\partial x_{t}} \approx {{\pi_{1}{\sum\limits_{i = {- 500}}^{500}{\left( {{i\; \Delta_{NOx}} + {y_{{MAP},{NOx}}\left( x_{t} \right)}} \right)\frac{e^{\frac{- {({i\; \Delta_{NOx}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}\Delta_{NOx}\frac{\partial{y_{{MAP},{NOx}}\left( x_{t} \right)}}{\partial x_{t}}}}} + {\frac{\pi_{2}}{2}{\sum\limits_{i = {- 500}}^{500}{e^{\frac{({{i\; \Delta_{CO}} + {y_{{MAP},{CO}}{(x_{t})}} - b_{co}})}{s}}\frac{e^{\frac{- {({i\; \Delta_{CO}})}^{2}}{2\; {\sigma_{{MAP},{CO}}^{2}{(x_{t})}}}}}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}\Delta_{CO}\frac{\partial{y_{{MAP},{CO}}\left( x_{t} \right)}}{\partial x_{t}}}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 136}\end{matrix}$

Using Equation 116 and Equation 117, it is found that:

$\begin{matrix}{\frac{\partial{J\left( x_{t} \right)}}{\partial x_{t}} \approx {{\pi_{1}{\sum\limits_{i = {- 500}}^{500}{\left( {{i\; \Delta_{NOx}} + {y_{{MAP},{NOx}}\left( x_{t} \right)}} \right)^{2}\frac{e^{\frac{- {({i\; \Delta_{NOx}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}\Delta_{NOx}\frac{\partial{{NN}\left( {x_{t},w_{{MAP},{NOx}}} \right)}}{\partial x_{t}}}}} + {\pi_{2}{\sum\limits_{i = {- 500}}^{500}{\frac{e^{\frac{({{i\; \Delta_{CO}} + {y_{{MAP},{CO}}{(x_{t})}} - b_{co}})}{s}}}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({i\; \Delta_{CO}})}^{2}}{2{\sigma_{{MAP},{CO}}^{2}{(x_{t})}}}}\Delta_{CO}\frac{\partial{{NN}\left( {x_{t},w_{{MAP},{NOx}}} \right)}}{\partial x_{t}}}}} + {\frac{\pi_{3}}{2}\left( {x_{t} - x_{t - 1}} \right)}}} & {{Equation}\mspace{14mu} 137}\end{matrix}$

It should be noticed that the form of the previous equation is verysimilar to the deterministic cost function of Equation 121. LikeEquation 121, the derivative of the first term in the equation above canbe computed by backpropagating the value:

$\left( {{i\; \Delta_{NOx}} + {y_{{MAP},{NOx}}\left( x_{t} \right)}} \right)\frac{\pi_{1}}{\sqrt{2\pi \; {\sigma_{{MAP},{NOx}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({i\; \Delta_{NOx}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}\Delta_{NOx}$

through the NOx neural network NN(x_(t),w_(MAP,NOx)) to the inputs.Likewise, the second term can be computed by backpropagating the value:

$\frac{\pi_{2}e^{\frac{({{i\; \Delta_{CO}} + {y_{{MAP},{CO}}{(x_{t})}} - b_{co}})}{s}}}{\sqrt{2\pi \; {\sigma_{{MAP},{CO}}^{2}\left( x_{t} \right)}}}e^{\frac{- {({i\; \Delta_{CO}})}^{2}}{2{\sigma_{{MAP},{NOx}}^{2}{(x_{t})}}}}\Delta_{CO}$

through the CO neural network NN(x_(t),w_(MAP,CO)) to the inputs. Thus,an efficient algorithm is available for the computation of both the costfunction and the derivative of the cost function. Using thesealgorithms, a gradient based approach with bounded values is used tosolve to the optimization problem with a probabilistic model. Theresulting optimized input vector can be output to the DCS andsubsequently used to control the plant to minimize NOx while maintainingCO below a threshold.

The above optimization problem shows one embodiment of the invention.However, there are many different approaches that could be used toinclude a probabilistic model in the optimization problem. In addition,the overall model confidence calculation and the confidence calculationper input could be included in an optimization problem. For example, inEquation 133, the overall model confidence of the NOx model could bemultiplied by the first term and the overall model confidence of the COmodel could be multiplied by the second term in the cost function. Thus,the model confidence could be used to adjust the cost function toemphasize components with more confident (accurate) models.

Optimization and Automatically Updated Models

The preferred embodiment of the invention is to use the trainingalgorithm in FIG. 8, limited to 3 cycles through the main loop, to trainthe neural network models of NOx and CO. Once the models are trained,they are used in an optimization system shown in FIG. 2 with theoptimization problem defined as given by Equation 131 and Equation 135.(The goal of the optimization problem is to reduce NOx, keep CO below aspecified bound, minimize the MV movement and maintain the MVs with aset of bounds. The optimization problem is solved every 15 minutes usinga NLP solver. The resulting optimized MVs are output to the DCS and heldconstant until the next optimization cycle.

In the preferred embodiment, the models are retrained once a day usingthe algorithm shown in FIG. 11. When data is collected, the optimizer isturned off, and the Design of Experiments with selected inputs is usedto collect the data. Once the additional data is collected, the model isretrained and the optimizer is turned back on using the updated models.It should be noted that if during data collection it is found that theMVs or CVs go beyond user specified bounds or the unit no longer is heldat a steady load, the data collection can be discontinued and theoptimizer may be re-enabled. Once the MVs or CVs are within bounds andthe unit is at a steady load, the data collection can be re-enabled andthe optimizer disabled until data collection is completed. Thus, thepreferred embodiment uses Bayesian Disturbance Rejection models that areupdated on-line in the optimization problem.

The preferred embodiment described above is used to control a set of MVsto reduce NOx emissions from a coal-fired power plant while maintainingCO below a limit. It should be noted that alternate embodiments of thisinvention exist for use in a power generation plant. For example, theinvention may be used to optimize the removal of soot in a boiler bycontrolling the time and execution of sootblowers. The BayesianDisturbance Rejection algorithm may be used to train models that predictthe effects of sootblower activations (MVs) on heating surfacecleanliness, temperature changes and heat duties. Given models of theeffect of sootblower activations, an optimization problem could becreated that is used to maximize the cleaning effects while observingconstraints on the blower activations. Thus, the invention could be usedfor sootblower optimization at a power generation unit.

As will be appreciated, the invention described herein may be useful inapplications beyond the power generation industry. For example, andwithout limitation, the invention may be used in the followingindustries: refining, chemicals, pulp and paper, pharmaceuticals, miningand manufacturing, and the food industries.

Yet another embodiment of the current invention includes its use fordynamic optimization rather than steady state optimization. Dynamicoptimization is described in detail in S. Piche, B. Sayyar-Rodsari, D.Johnson and M. Gerules, “Nonlinear model predictive control using neuralnetworks,” IEEE Control Systems Magazine, vol. 20, no. 2, pp. 53-62,2000, which is fully incorporated herein by reference. As will beappreciated, dynamic optimization relies on the use of a dynamic model.In Equation 1, the model is defined as:

m _(t)=NN(x _(t−1) ,w)  Equation 138

Alternatively, the model could be defined as:

m _(t)=NN(x _(t−1) ,w ¹)+NN(x _(t−2) ,w ²)+ . . . +NN(x _(t−M) ,w^(M))  Equation 139

where M is the number of neural network models. If the sample period ofthe model of the equation above is significantly less than the steadystate response period, typically 1/10^(th) of the steady state response,this type of model forms a finite impulse response model and can be usedto predict dynamics of the process. By extension, a Bayesian DisturbanceRejection algorithm can be used for training the M models of theprevious equation. Once trained, these models could be used in amulti-step ahead optimization scheme. Thus, the neural model could beused for trajectory based model optimization. For example, models couldbe used to predict over the horizon from time (t+1) to t+M and then theinputs could be optimized over the horizon from (t) to t+M−1 given theappropriate cost function and constraints. In the preferred embodiment,the case was shown where M=1, but all results and conclusions providedhere can be extended to cases where M>1.

A neural network based model is used in the preferred embodiment form_(t). According to other possible embodiments, other empirical modelforms may be used including, for example, a linear model, polynomialcurve fit, a least squares support vector machine, a support vectormachine, and kernel methods including Gaussian Processes.

Thus, exemplary embodiments of the present invention include a methodand/or system for training a disturbance rejection model that isconfigured to model an operation of a system so to calculate a predictedvalue for a system output at a future time. The disturbance rejectionmodel may include a network for mapping system inputs to the systemoutput, the network including a weight vector and a feedbackcoefficient. The method may include the steps of: obtaining a trainingdataset, the training dataset including a time series dataset of actualvalues for the system inputs and the system output as determined frommeasurements taken of the operation of the system during a historicaloperating period and training the disturbance rejection model pursuantto the training dataset. The training may include calculating updatedvalues for each of the weight vector and the feedback coefficient of thenetwork by minimizing an error function that include a firsthyperparameter and a second hyperparameter. The first hyperparameter mayinclude a vector for penalizing the weight vector and the secondhyperparameter may include a scalar.

As described, the disturbance rejection model may include a disturbancerejection configuration in which the predicted value made by thedisturbance rejection model for the system output at the future time isbased upon: a predicted value made by the network for the system outputat the future time and a value of a bias, the bias being based upon thefeedback coefficient and an error. The bias may include the errormultiplied by the feedback coefficient, The error includes a differencebetween: a predicted value made by the neural network of the systemoutput at the previous time and an actual value of the system output atthe previous time, The actual value is based upon a measurement taken bya sensor disposed in the system for measuring an operating parameterthat relates to the system output.

Further, the network may be a neural network that includes multiplelayers having nodes, the multiple layers including at least an inputlayer, an output layer, one or more hidden layers, and forward weightmatrixes. The input layer may include a plurality of the nodes, theplurality of the nodes corresponding respectively to the system inputs.Each of the plurality of the nodes may be configured to receive an inputsignal relating to a value for a particular system input. The outputlayer may include at least one of the nodes, which corresponds to thesystem output. The one or more hidden layers may be disposed between theinput layer and the output layer, each of the one or more hidden layersincluding a plurality of the nodes. The forward weight matrices mayinclude connectors that connect the nodes of successive layers of themultiple layers of the neural network and a weight value for each of theconnectors. The weight vector may define the weight values for theconnectors of the forward weight matrices.

The training dataset may be generated via a design of experimentprocess. The design of experiment process may include the manipulatedvariables being moved independently such that moves resulting therefromare uncorrelated over time. The training dataset may include T samples,the T samples taken at times=[1, 2, . . . , t, . . . , T], wherein time(t), time (t−1), time (t−2) denote generalized samples within thetraining dataset, where the time (t−2) occurs just prior to the time(t−1) and the time (t−1) occurs just prior to the time (t).

The system be a power generating unit that includes at least a steamturbine operably connected to a boiler. The system output may include acontrolled variable, and the system inputs may include manipulatedvariables and disturbance variables. The manipulated variables mayinclude any of the following: an over-fired air damper position, a tiltbias, a yaw bias, a secondary air damper position bias, a fuel millspeed bias, a primary air bias, an O2 bias, and a master burner tilt.The controlled variable may include an emission level for the powergenerating unit.

The training may include: a first step that is defined as thecalculating of the updated values for the weight vector and the feedbackcoefficient by minimizing the error function; and a second step that mayinclude calculating updated values for the first hyperparameter and thesecond hyperparameter based upon the updated values that are calculatedin the first step for the weight vector and the feedback coefficient.The computing the updated values for each of the weight vector and thefeedback coefficient may include a matrix version of the disturbancerejection model. The step of computing the updated values for the weightvector and the feedback coefficient that minimize the error function mayinclude a first iterative process, each iteration thereof beinginclusive of the following steps: using forward propagation to computethe vector y using the matrix version of the disturbance rejectionmodel; using backward propagation on the disturbance rejection model tocompute a derivative of the error function with respect to the weightvector and the feedback coefficient; and using a gradient basedalgorithm to compute a current iteration value for the weight vectorwhile constraining a current iteration value of the feedback coefficientto values between 0.0 and 1.0. The gradient based algorithm may includea steepest descent algorithm.

The first iterative process for computing the updated values for theweight vector and the feedback coefficient may continue until anultimate iteration is found that results in a change in the value of theerror function from a preceding iteration to be less than a predefinedthreshold. The current iteration value for each of the weight vector andthe feedback coefficient found in the ultimate iteration may be deemedto be the values of the weight vector and the feedback coefficient thatminimize the error function and, thus, are designated as the updatedvalues for the weight vector and the feedback coefficient.

The step of calculating the re-estimated value for the firsthyperparameter and the second hyperparameter may include a seconditerative process. A single iteration of the second iterative processmay include the steps of: calculating of the updated values for theweight vector and the feedback coefficient given the re-estimated valuesof the first hyperparameter and the second hyperparameter as calculatedin an iteration of the second iterative process that is previous to thesingle iteration; and calculating the re-estimated value for each of thefirst hyperparameter and the second hyperparameter given the updatedvalues for the weight vector and the feedback coefficient as calculatedwithin the single iteration. The second iterative process of thedisturbance rejection training algorithm may continue until there-estimated value for each of the first hyperparameter and the secondhyperparameter converges to a value.

The present invention further describes a related system that includes apower generating unit and a control system operably connected to thepower generating unit for controlling an operation thereof. The controlsystem may include: a hardware processor; and a machine readable storagemedium on which is stored the disturbance rejection model andinstructions that cause the hardware processor to execute a processrelated to control of the power generating unit. The process mayinclude: obtaining a training dataset; and training the disturbancerejection model pursuant to the training dataset. The training mayinclude calculating updated values for each of the weight vector and thefeedback coefficient of the neural network by minimizing an errorfunction that includes a first hyperparameter and a secondhyperparameter. The first hyperparameter may include a vector forpenalizing the weight vector and the second hyperparameter may include ascalar.

According to another related exemplary embodiment, the present inventionfurther describes a method for controlling an operation of a system thatincludes a disturbance rejection model that is configured for modelingthe operation of the system so to generate a predicted value for asystem output at a future time. The disturbance rejection model mayinclude a neural network for mapping system inputs to the system output.The method may include the steps of: calculating a probabilisticdistribution for the predicted value of the system output at the futuretime, where the calculating of the probabilistic distribution includes aBayesian evidence framework without sampling; and controlling theoperation of the system per the calculated probabilistic distribution.

The calculating of the probabilistic distribution may includecalculating a mean (μ) and a variance (σ̂2) of a normal distribution (N).The mean (μ) may correspond to the predicted value made by the traineddisturbance rejection model for the system output at the future time.The variance (σ̂2) may depend, at least in part, on a value of an inverseof the second hyperparameter.

The controlling the operation of the system may include optimizing theoperation of the system. The system may include a power generating unitthat includes at least a steam turbine operably connected to a boiler.The manipulated variables may include any of the following: anover-fired air damper position, a tilt bias, a yaw bias, a secondary airdamper position bias, a fuel mill speed bias, a primary air bias, an O2bias, a master burner tilt, a corner burner bias, and a windbox tofurnace differential pressure bias. The controlled variable may include:an emission level for the power generating unit; a heat rate for thepower generating unit; or a boiler efficiency for the power generatingunit. Alternatively, the controlled variable may include at least oneof: a reheat temperature level for the power generating unit; asuperheat steam temperature level for the power generating unit; areheat spray level for the power generating unit; and a superheat steamspray level for the power generating unit. The boiler may include a heattransfer device. In such cases, the manipulated variables may include acumulative time since a last activation of a first or second sootblower. The controlled variable may also include one of the followingaspects associated with the heat transfer device: a cleanliness factor;a heat transfer coefficient; and a heat duty. The controlled variablemay also include a steam inlet and output temperature differentialassociated with a heat transfer device in the boiler of the powergenerating unit.

The optimizing the operation of the system may include the steps of:determining one or more operating goals associated with the operation ofthe power generating unit; determining one or more operating constraintsassociated with the operation of the power generating unit; andproviding the one or more operating goals and the one or more operatingconstraints to an optimization system. The optimization system mayinclude: the disturbance rejection model; an optimizer for determiningoptimal setpoint values for the manipulated variables associated withthe control of the power generating unit, the optimal setpoint valuesdetermined in accordance with the one or more operating goals and theone or more operating constraints; and a communication link to a controlsystem for communicating the optimal setpoint values thereto. Thecontrol system may be configured to provide closed-loop control for theoperation of the power generating unit. The optimizing may furtherinclude the steps of: using the optimization system and theprobabilistic distribution for the predicted value of the system outputat the future time as calculated by the disturbance rejection model tocalculate the optimal setpoint values for operating the power generatingunit; and communicating via the communication link the optimal setpointvalues to the control system of the power generating unit. The optimizerof the optimization system may be selected from a group consisting of:linear programming, quadratic programming, mixed integer non-linearprogramming, stochastic programming, global non-linear programming,genetic algorithms, and particle/swarm techniques.

The method may further include the step of receiving a cost functionthat defines costs associated with each of the one or more operatinggoals and the one or more operating constraints. The optimizer of theoptimization system may be configured to determine the optimal setpointvalues for the manipulated variables by accessing the disturbancerejection model to minimize the cost function. The method may furtherinclude instructing the control system to control the power generatingunit in accordance with the optimal setpoints once the optimal setpointsare communicated thereto. At least one of the costs defined within thecost function may be configured to vary according to a value of anoperating goal probability. The operating goal probability may define aprobability as to whether a selected one of the one or more operatinggoals will be satisfied at the future time. The goal probability may bederived from the probabilistic distribution for the predicted value ofthe system output at the future time. The operating goals may includeany of the following: an output level for the power generating unit; anda heat rate for the power generating unit. At least one of the costsdefined within the cost function may be configured to vary according toa value of an operating constraint probability. The operating constraintprobability may define a probability as to whether a selected one of theone or more operating constraints is violated. The operating constraintprobability may be derived from the probabilistic distribution for thepredicted value of the system output at the future time. The operatingconstraints may include any of the following: a limit on emissions ofNOx a limit on emissions of SO2 a limit on emissions of CO2 and a limiton emissions of CO. The cost function may include a mathematicalrepresentation for evaluating the operation of the power generating unitrelative to the one or more operating goals and the one or moreoperating constraints. The expression of the one or more operating goalsin the cost function may include coefficients for establishing arelative weighting therebetween.

Another exemplary embodiment of the present invention includes a methodfor modeling an operation of a system that may include a disturbancerejection model that is configured to generate a predicted value for asystem output at a future time. The disturbance rejection model mayinclude a neural network for mapping system inputs to the system output.The method may include the steps of: training the disturbance rejectionmodel per a training dataset; and calculating a confidence metric forthe disturbance rejection model. As described herein, the confidencemetric is configured to indicate a probability that a predicted sign ofa gain in the system output at the future time made by the disturbancerejection model is correct. The confidence metric may be based on theprobabilistic distribution discussed above.

The algorithm for calculating the confidence metric may be configuredsuch that the confidence metric may include: a high probability that thepredicted sign of the gain in the system output at the time (t) iscorrect when the variance of the probabilistic distribution is smallcompared to a difference between: the predicted value of the systemoutput made by the trained disturbance rejection model at the time (t)and the system output (d_(t−1)) at the time (t−1); and a low probabilitythat the predicted sign of the gain in the system output at the time (t)is correct when the variance of the probabilistic distribution is largecompared to the difference between: the predicted value of the systemoutput made by the trained disturbance rejection model at the time (t)and the system output (d_(t−1)) at the time (t−1).

The predicted sign of the gain in the system output is defined as apredicted sign of a value change for the system output between the time(t−1) and the time (t). The predicted sign may include: a positive valuewhen the value change results in the predicted value for the systemoutput at the time (t) increasing from the actual value of the systemoutput (d_(t−1)) at the time (t−1); and a negative value when the valuechange results in the predicted value for the system output at the time(t) decreasing from the actual value of the system output (d_(t−1)) atthe time (t−1).

The calculation of a confidence score based on the confidence metric forthe disturbance rejection model may include: calculating the confidencemetric for selected samples chosen from the T samples of the trainingdataset; calculating an average value for the confidence metric for theselected samples; and designating the calculated average value as theconfidence score for the disturbance rejection model. The confidencescore may include a percentage indicating a likelihood that thepredicted sign of the gain in the system output at the future time madeby the disturbance rejection model is correct based upon the performanceof a whole of the disturbance rejection model over the selected samplesfrom the training dataset. The future time may include a time that issubsequent to the historical operating period on which the trainingdataset is based. The confidence score may include a spectrum withinwhich: a value of 100% indicates that the disturbance rejection modelpredicts the correct sign of the gain for 100% of the samples within theselected samples from the training dataset; and a value of 50% indicatesthat the disturbance rejection model has no confidence in predicting thecorrect sign of the gain for the selected samples from the trainingdataset.

The method may further include the step of recommending a modificationto the disturbance rejection model based upon a comparison of theconfidence score to a predetermined minimum threshold. The modificationmay relate to a need for augmenting the training dataset of thedisturbance rejection model with additional data when the confidencescore fails to satisfy the predetermined minimum threshold.

A related exemplary embodiment includes a method for modeling a system,the method including a disturbance rejection model configured formodeling an operation of the system so to generate a predicted value fora system output at a future time. The disturbance rejection model mayinclude a neural network for mapping system inputs to the system outputand input-output pairings, each of which represent a unique pairing ofone of the system inputs with the system output. The method may includethe steps of: calculating a confidence metric for a selected one of theinput-output pairings of the disturbance rejection model; andrecommending a modification be made to the disturbance rejection modelbased upon the confidence metric calculated for the selected one of theinput-output pairing. The confidence metric may be configured toindicate a probability that a predicted sign of a gain in the systemoutput at the future time made by the disturbance rejection model iscorrect when the system input of the selected one of the input-outputpairings is varied.

In accordance with the present example, adjacent samples within the Tsamples of the training dataset are defined as paired samples thatcorrespond to consecutively occurring times within the times=[1, 2, . .. , t, . . . , T]. Pertaining to the training dataset, a move is definedas a change to a system input setting that occurs between theconsecutively occurring times of one of the adjacent samples. Thetraining dataset may be configured such that each of the adjacentsamples defined within the training dataset include no more than one ofthe moves. Alternatively, the training dataset may be configured suchthat most of the adjacent samples defined within the training datasetinclude more than one of the moves.

The confidence metric may be used to calculate a confidence score forthe selected one of the input-output pairings. The calculation of theconfidence score may include: determining a subset of samples from the Tsamples of the training dataset, the subset of samples may include onlythose adjacent samples in which the move involves the selected one ofthe input-output pairings; calculating the confidence metric for each ofthe samples within the subset of samples; calculating an average valuefor the confidence metric for the subset of samples; and designating theaverage value as the confidence score for the selected one of theinput-output pairings. The confidence score determined for the selectedone of the input-output pairings may be compared to a predeterminedthreshold. The recommended modification for the disturbance rejectionmodel may be based on the comparison of the confidence score to thepredetermined threshold. The recommended modification may relate towhether additional data should be obtained for the training dataset. Thepredetermined threshold may include a minimum threshold. The recommendedmodification for the disturbance rejection model may include at leastone of the following: when the comparison shows that the confidencescore fails to satisfy the minimum threshold, the recommendedmodification may include one regarding a high need for obtaining theadditional data for training dataset; and when the comparison shows thatthe confidence score satisfies the minimum threshold, the recommendedmodification may include one regarding a low need for obtaining theadditional data for training dataset. The recommended modification forthe disturbance rejection model may include defining the additional dataas data corresponding to the system input of the selected one of theinput-output pairings.

The calculation of the confidence score may be repeated for at least oneother of the input-output pairings such that the confidence score isdetermined for two or more of the input-output pairings. The confidencescores determined for the two or more of the input-output pairings areeach compared against a predetermined threshold. The recommendedmodification for the disturbance rejection model may include one relatedto a need to obtain additional data for the training dataset based uponhow each of the confidence scores compared to the predeterminedthreshold. Alternatively, the confidence scores determined for the twoor more of the input-output pairings may be compared against each otherso to determine at least a first input-output pairing that has aconfidence score that is preferable compared to a second input-outputpairing. The recommended modification for the disturbance rejectionmodel may include one in which the need to obtain the additional datafor the system input of the second input-output pairing is prioritizedover the need to gather additional data for the system input of thefirst input-output pairing.

As one of ordinary skill in the art will appreciate, the many varyingfeatures and configurations described above in relation to the severalexemplary embodiments may be further selectively applied to form theother possible embodiments of the present invention. For the sake ofbrevity and taking into account the abilities of one of ordinary skillin the art, all of the possible iterations is not provided or discussedin detail, though all combinations and possible embodiments embraced bythe several claims below or otherwise are intended to be part of theinstant application. In addition, from the above description of severalexemplary embodiments of the invention, those skilled in the art willperceive improvements, changes and modifications. Such improvements,changes and modifications within the skill of the art are also intendedto be covered by the appended claims. Further, it should be apparentthat the foregoing relates only to the described embodiments of thepresent application and that numerous changes and modifications may bemade herein without departing from the spirit and scope of theapplication as defined by the following claims and the equivalentsthereof.

That which is claimed:
 1. A method for controlling an operation of asystem that includes a disturbance rejection model that is configuredfor modeling the operation of the system so to generate a predictedvalue for a system output at a future time, the disturbance rejectionmodel comprising a network for mapping system inputs to the systemoutput, the method comprising the steps of: calculating a probabilisticdistribution for the predicted value of the system output at the futuretime, wherein the calculating of the probabilistic distributioncomprises a Bayesian evidence framework without sampling; andcontrolling the operation of the system per the calculated probabilisticdistribution for the predicted value for the system output at the futuretime.
 2. The method of claim 1, further comprising the step of trainingthe disturbance rejection model according to a disturbance rejectiontraining algorithm and a training dataset; wherein the training datasetcomprising a time series dataset of actual values for the system inputsand the system output as determined from measurements taken of theoperation of the system during a historical operating period, thehistorical operating period occurring prior to the future time; andwherein the disturbance rejection model comprises a disturbancerejection configuration in which the predicted value made by thedisturbance rejection model for the system output at the future time isbased upon: a predicted value made by the network for the system outputat the future time; and a value of a bias that is based upon a feedbackcoefficient and an error.
 3. The method of claim 2, wherein the networkcomprises a neural network that includes multiple layers having nodes,the multiple layers including at least an input layer, an output layer,one or more hidden layers, and forward weight matrixes; and wherein: theinput layer comprises a plurality of the nodes, the plurality of thenodes corresponding respectively to the system inputs, wherein each ofthe plurality of the nodes is configured to receive an input signalrelating to a value for a particular one of the system inputs; theoutput layer comprises at least one of the nodes, the at least one ofthe nodes corresponding to the system output; the one or more hiddenlayers are disposed between the input layer and the output layer, eachof the one or more hidden layers comprising a plurality of the nodes;and the forward weight matrices comprise connectors that connect thenodes of successive layers of the multiple layers of the neural networkand a weight value for each of the connectors; wherein a weight vectordefines the weight values for the connectors of the forward weightmatrices; and wherein the bias comprises the error multiplied by thefeedback coefficient, and wherein the error comprises a differencebetween: a predicted value made by the neural network of the systemoutput at the previous time; and an actual value of the system output atthe previous time, wherein the actual value is based upon a measurementtaken by a sensor disposed in the system for measuring an operatingparameter related to the system output.
 4. The method of claim 3,wherein the calculating of the probabilistic distribution comprisescalculating a mean (μ) and a variance (σ²) of a normal distribution (

), wherein: the mean (μ) corresponds to the predicted value made by thetrained disturbance rejection model for the system output at the futuretime; and the variance (σ²) depends, at least in part, on a value of aninverse of the second hyperparameter.
 5. The method of claim 3, wherein:the future time is denoted as a future time (t); the previous time,which is defined as just preceding the future time (t), is denoted as afirst previous time (t−1); and a second previous time, which is definedas just preceding the first previous time (t−1), is denoted as a secondprevious time (t−2); wherein the system inputs at the first previoustime (t−1) and the second previous time (t−2) are described by inputvectors, which are denoted, respectively, as a first input vector(x_(t−1)) and a second input vector (x_(t−2)); wherein the system outputat the first previous time (t−1) is denoted as a system output(d_(t−1)); wherein the step of calculating the probabilisticdistribution for the predicted value of the system output at the futuretime comprises calculating a mean (μ) and a variance (σ²) that define anormal distribution (

); and wherein each of the first previous time (t−1) and the secondprevious time (t−2) comprises times that are subsequent to thehistorical operating period on which the training dataset is based. 6.The method of claim 5, wherein the mean (μ) comprises the predictedvalue made by the disturbance rejection model for the system output atthe future time (t), as given by:μ(x _(t−1) ,x _(t−2) ,d _(t−1) ,k*,w _(MAP))=y(x _(t−1) ,x _(t−2) ,d_(t−1) ,k*,w _(MAP)) where: the y_(t) comprises the predicted value madeby the disturbance rejection model for the system output at the futuretime (t); the k* denotes a trained value for the feedback coefficient(k); and the w_(MAP) denotes a trained value of the weight vector;wherein each of the trained values, the k* and the w_(MAP), aredetermined from the disturbance rejection model once trained via thedisturbance rejection training algorithm pursuant to the trainingdataset.
 7. The method of claim 6, wherein the variance (σ²) is definedusing the Bayesian evidence framework as:σ²(x _(t−1) ,x _(t−2) ,d _(t−1) ,k*,w _(MAP))=β⁻¹ +g _(t) ^(T) A ⁻¹ g_(t) wherein: the A is a second derivative of a negative log posteriordistribution as given by:A=Diag(α)I+βH; where: the I comprises an identity matrix; and the Hcomprises a Hessian matrix, the Hessian matrix being computed using aR-propagation algorithm; and the g_(t) comprises a vector g_(t) that iscalculated using at least one of the following equations:${{{{{g_{t} = \frac{\partial{y\left( {x_{t - 1},x_{t - 2},k,w} \right)}}{\partial w}}}_{{k = k^{*}},{w = w_{MAP}}};{and}}{g_{t} = \begin{pmatrix}{\frac{\partial{{NN}\left( {x_{t - 1},w} \right)}}{\partial w} - k} & \frac{\partial{{NN}\left( {x_{t - 2},w} \right)}}{\partial w}\end{pmatrix}}}}_{{k = k^{*}},{w = w_{MAP}}}.$
 8. The method of claim6, wherein the system output comprises a controlled variable; andwherein the system inputs comprise disturbance variables and manipulatedvariables, wherein: the manipulated variables are defined as variablesthat a system operator is able to manipulate to affect the controlledvariable; and the disturbance variables are defined as variables thataffect the controlled variable that the system operator is unable tomanipulate; wherein the controlling the operation of the system includesoptimizing the operation of the system.
 9. The method of claim 8,wherein the system comprises a power generating unit, the powergenerating unit comprising at least a steam turbine operably connectedto a boiler; and wherein the manipulated variables comprise two or moreof the following: an over-fired air damper position, a tilt bias, a yawbias, a secondary air damper position bias, a fuel mill speed bias, aprimary air bias, an O₂ bias, a master burner tilt, a corner burnerbias, and a windbox to furnace differential pressure bias.
 10. Themethod of claim 9, wherein the controlled variable comprises at leastone of: an emission level for the power generating unit; a heat rate forthe power generating unit; and a boiler efficiency for the powergenerating unit.
 11. The method of claim 9, wherein the controlledvariable comprises at least one of: a reheat temperature level for thepower generating unit; and a superheat steam temperature level for thepower generating unit; a reheat spray level for the power generatingunit; and a superheat steam spray level for the power generating unit.12. The method of claim 8, wherein the system comprises a powergenerating unit, the power generating unit comprising at least a steamturbine operably connected to a boiler, wherein the boiler comprises aheat transfer device; wherein the manipulated variables comprise atleast one of: a cumulative time since a last activation of a first sootblower; and a cumulative time since the last activation of a second sootblower; wherein the controlled variable comprises at least one of thefollowing aspects associated with the heat transfer device: acleanliness factor; a heat transfer coefficient; and a heat duty. 13.The method of claim 8, wherein the system comprises a power generatingunit, the power generating unit comprising at least a steam turbineoperably connected to a boiler; wherein the manipulated variablescomprise at least one of: a cumulative time since a last activation of afirst soot blower; and a cumulative time since the last activation of asecond soot blower; wherein the controlled variable comprises a steaminlet to steam outlet differential temperature associated with a heattransfer device in the boiler of the power generating unit.
 14. Themethod of claim 9, wherein the optimizing the operation of the systemcomprises the steps of: determining one or more operating goalsassociated with the operation of the power generating unit; determiningone or more operating constraints associated with the operation of thepower generating unit; providing the one or more operating goals and theone or more operating constraints to an optimization system thatincludes: the disturbance rejection model; an optimizer for determiningoptimal setpoint values for the manipulated variables associated withthe control of the power generating unit, the optimal setpoint valuesdetermined in accordance with the one or more operating goals and theone or more operating constraints; and a communication link to a controlsystem for communicating the optimal setpoint values thereto, whereinthe control system is configured to provide closed-loop control for theoperation of the power generating unit; using the optimization systemand the probabilistic distribution for the predicted value of the systemoutput at the future time as calculated by the disturbance rejectionmodel to calculate the optimal setpoint values for operating the powergenerating unit; and communicating via the communication link theoptimal setpoint values to the control system of the power generatingunit.
 15. The method of claim 14, wherein the optimizer of theoptimization system is selected from a group consisting of: linearprogramming, quadratic programming, mixed integer non-linearprogramming, stochastic programming, global non-linear programming,genetic algorithms, and particle/swarm techniques; further comprisingthe steps of: receiving a cost function that defines costs associatedwith each of the one or more operating goals and the one or moreoperating constraints, wherein the optimizer of the optimization systemis configured to determine the optimal setpoint values for themanipulated variables by accessing the disturbance rejection model tominimize the cost function; and instructing the control system tocontrol the power generating unit in accordance with the optimalsetpoints once the optimal setpoints are communicated thereto.
 16. Themethod of claim 15, wherein at least one of the costs defined within thecost function is configured to vary according to a value of an operatinggoal probability, the operating goal probability defining a probabilityas to whether a selected one of the one or more operating goals will besatisfied at the future time; and wherein the goal probability isderived from the probabilistic distribution for the predicted value ofthe system output at the future time; wherein the selected one of theone or more operating goals comprises at least one of: an output levelfor the power generating unit; and a heat rate for the power generatingunit.
 17. The method of claim 15, wherein at least one of the costsdefined within the cost function is configured to vary according to avalue of an operating constraint probability, the operating constraintprobability defining a probability as to whether a selected one of theone or more operating constraints is violated; wherein the operatingconstraint probability is derived from the probabilistic distributionfor the predicted value of the system output at the future time; andwherein the selected one of the one or more operating constraintscomprises one of the following: a limit on emissions of NOx; a limit onemissions of SO₂; a limit on emissions of CO₂; and a limit on emissionsof CO.
 18. The method of claim 15, wherein the cost function comprises amathematical representation for evaluating the operation of the powergenerating unit relative to the one or more operating goals and the oneor more operating constraints; and wherein an expression of the one ormore operating goals in the cost function include coefficients forestablishing a relative weighting therebetween.
 19. The method of claim9, wherein the disturbance rejection training algorithm comprisescalculating updated values for each of the weight vector and thefeedback coefficient of the network by minimizing an error function thatincludes a first hyperparameter and a second hyperparameter; and whereinthe first hyperparameter includes a vector for penalizing the weightvector and the second hyperparameter comprises a scalar.
 20. The methodof claim 19, wherein the training dataset comprises T samples, the Tsamples taken at times=[1, 2, . . . , t, . . . , T], wherein a time (t),a time (t−1), a time (t−2) denoted generalized samples in the trainingdataset where the time (t−2) occurs just prior to the time (t−1) and thetime (t−1) occurs just prior to the time (t); wherein the computing theupdated values for each of the weight vector and the feedbackcoefficient includes a matrix version of the disturbance rejection modelin which a C matrix comprises a T by T matrix configured to represent acoupling between time series elements of the training dataset; whereinthe error function, which is denoted as E below, that is minimized tocalculate the updated values of the weight vector and the feedbackcoefficient is defined as:E=½w ^(T)(Diag(α))w+½βε^(T) Cε wherein: the C comprises the C matrix;the α comprises the first hyperparameter; the β comprises the secondhyperparameter; the w^(T) comprises a transposition of the weight vector(w); the Diag(α) comprises a matrix with a diagonal element equal to avalue of a corresponding diagonal element of the vector of the firsthyperparameter while all other values in the matrix of the Diag(α) areset to a value of 0; the ε comprises an error vector that comprises adifference between the predicted values of the system output by thedisturbance rejection model and the corresponding actual values of thesystem output at the times of the training dataset; and the ε^(T)comprises a transposition of the error vector.
 21. A system comprising:a power generating unit; a control system operably connected to thepower generating unit for controlling an operation thereof, the controlsystem comprising: a hardware processor; and a machine readable storagemedium on which is stored: a disturbance rejection model; andinstructions that cause the hardware processor to execute a processrelated to control of the power generating unit; wherein the disturbancerejection model models the operation of the power generating unit so togenerate a predicted value for an output of the power generating unit ata future time, the disturbance rejection model comprising: a network formapping inputs of the power generating unit to the output; wherein theprocess comprises: calculating a probabilistic distribution for thepredicted value of the output at the future time, wherein thecalculating of the probabilistic distribution comprising a Bayesianevidence framework without sampling; and controlling the operation ofthe power generating until per the calculated probabilistic distributionfor the predicted value for the output at the future time.
 22. Thesystem of claim 21, wherein the process comprises the step of trainingthe disturbance rejection model according to a disturbance rejectiontraining algorithm and a training dataset; wherein the training datasetcomprising a time series dataset of actual values for the inputs and theoutput as determined from measurements taken of the operation of thepower generating unit during a historical operating period, thehistorical operating period occurring prior to the future time; andwherein the disturbance rejection model comprises a disturbancerejection configuration in which the predicted value made by thedisturbance rejection model for the output at the future time is basedupon: a predicted value made by the network for the output at the futuretime; and a value of a bias that is based upon a feedback coefficientand an error.
 23. The system of claim 22, wherein the network comprisesa neural network that includes multiple layers having nodes, themultiple layers including at least an input layer, an output layer, oneor more hidden layers, and forward weight matrixes; and wherein: theinput layer comprises a plurality of the nodes, the plurality of thenodes corresponding respectively to the inputs, wherein each of theplurality of the nodes is configured to receive an input signal relatingto a value for a particular one of the inputs; the output layercomprises at least one of the nodes, the at least one of the nodescorresponding to the output; the one or more hidden layers are disposedbetween the input layer and the output layer, each of the one or morehidden layers comprising a plurality of the nodes; and the forwardweight matrices comprise connectors that connect the nodes of successivelayers of the multiple layers of the neural network and a weight valuefor each of the connectors; wherein a weight vector defines the weightvalues for the connectors of the forward weight matrices; and whereinthe bias comprises the error multiplied by the feedback coefficient, andwherein the error comprises a difference between: a predicted value madeby the neural network of the output at the previous time; and an actualvalue of the output at the previous time, wherein the actual value isbased upon a measurement taken by a sensor disposed in the powergenerating unit for measuring an operating parameter related to theoutput.
 24. The system of claim 23, wherein: the future time is denotedas a future time (t); the previous time, which is defined as justpreceding the future time (t), is denoted as a first previous time(t−1); and a second previous time, which is defined as just precedingthe first previous time (t−1), is denoted as a second previous time(t−2); wherein the inputs at the first previous time (t−1) and thesecond previous time (t−2) are described by input vectors, which aredenoted, respectively, as a first input vector (x_(t−1)) and a secondinput vector (x_(t−2)); wherein the output at the first previous time(t−1) is denoted as a output (d_(t−1)); wherein the step of calculatingthe probabilistic distribution for the predicted value of the output atthe future time comprises calculating a mean (μ) and a variance (σ²)that define a normal distribution (

); and wherein each of the first previous time (t−1) and the secondprevious time (t−2) comprises times that are subsequent to thehistorical operating period on which the training dataset is based. 25.The system of claim 24, wherein the mean (μ) comprises the predictedvalue made by the disturbance rejection model for the output at thefuture time (t), as given by:μ(x _(t−1) ,x _(t−2) ,d _(t−1) ,k*,w _(MAP))=y(x _(t−1) ,x _(t−2) ,d_(t−1) ,k*,w _(MAP)) where: the y_(t) comprises the predicted value madeby the disturbance rejection model for the output at the future time(t); the k* denotes a trained value for the feedback coefficient (k);and the w_(MAP) denotes a trained value of the weight vector; whereineach of the trained values, the k* and the w_(MAP), are determined fromthe disturbance rejection model once trained via the disturbancerejection training algorithm pursuant to the training dataset.
 26. Thesystem of claim 25, wherein the output comprises a controlled variableand the inputs comprise disturbance variables and manipulated variables;wherein the power generating unit comprises at least a steam turbineoperably connected to a boiler; wherein the controlling the operation ofthe power generating unit includes optimizing the operation of the powergenerating unit; wherein the system further comprises an optimizationsystem that includes: the disturbance rejection model; and an optimizerfor determining optimal setpoint values for the manipulated variablesassociated with the control of the power generating unit, the optimalsetpoint values determined in accordance with the one or more operatinggoals and the one or more operating constraints; wherein the processfurther includes: determining one or more operating goals associatedwith the operation of the power generating unit; determining one or moreoperating constraints associated with the operation of the powergenerating unit; providing the one or more operating goals and the oneor more operating constraints to the optimization system; using theoptimization system and the probabilistic distribution for the predictedvalue of the system output at the future time as calculated by thedisturbance rejection model to calculate the optimal setpoint values foroperating the power generating unit; and communicating the optimalsetpoint values to the control system of the power generating unit. 27.The system of claim 26, wherein the optimizer of the optimization systemis selected from a group consisting of: linear programming, quadraticprogramming, mixed integer non-linear programming, stochasticprogramming, global non-linear programming, genetic algorithms, andparticle/swarm techniques; wherein the process further comprises:receiving a cost function that defines costs associated with each of theone or more operating goals and the one or more operating constraints,wherein the optimizer of the optimization system is configured todetermine the optimal setpoint values for the manipulated variables byaccessing the disturbance rejection model to minimize the cost function;and instructing the control system to control the power generating unitin accordance with the optimal setpoints once the optimal setpoints arecommunicated thereto.
 28. The system of claim 27, wherein at least oneof the costs defined within the cost function is configured to varyaccording to a value of an operating goal probability, the operatinggoal probability defining a probability as to whether a selected one ofthe one or more operating goals will be satisfied at the future time;and wherein the goal probability is derived from the probabilisticdistribution for the predicted value of the output at the future time;wherein the selected one of the one or more operating goals comprises atleast one of: an output level for the power generating unit; and a heatrate for the power generating unit.
 29. The system of claim 27, whereinat least one of the costs defined within the cost function is configuredto vary according to a value of an operating constraint probability, theoperating constraint probability defining a probability as to whether aselected one of the one or more operating constraints is violated;wherein the operating constraint probability is derived from theprobabilistic distribution for the predicted value of the output at thefuture time; and wherein the selected one of the one or more operatingconstraints comprises one of the following: a limit on emissions of NOx;a limit on emissions of SO₂; a limit on emissions of CO₂; and a limit onemissions of CO.